Opened 4 years ago

Closed 4 years ago

#272 closed defect - wrong answer (fixed)

line shielding function is poorly behaved for large damping coefficient

Reported by: Gary J. Ferland Owned by: nobody
Priority: critical Milestone:
Component: radiative transfer Version: trunk
Keywords: Cc:

Description

We now use the Federman+79 line self shielding function for all lines in the code. It was derived for small damping coefficients and the damping wings are now treated properly for large damping coefficients.

http://cdsads.u-strasbg.fr/abs/1979ApJ...227..466F

used at rt_continuum_shield_fcn.cpp:63 see google groups thread under "convergence branch" dates 2013 nov 24 starting with Ryan email.

The relevant part of Federman+79 is the bit around A3 and A9 - they do not use standard stellar atmosphere (ie, Chandra) notation but gamma/beta is the ratio of damping to doppler widths, so r is our damping coefficient.

A9 handles the damping wings - looks like the function becomes large when r is large. From Fig 5 of Ferland+13 (Fe UTA) we can get a as large as 1000. we are likely saved by the min at line 111 - otherwise the shield would become quite large at small tau, large a.

Change History (3)

comment:1 Changed 4 years ago by Gary J. Ferland

ran auto at r8626 with attached patch, which caps a at 1e-2 - federman says approximation optimized for 1 = 1e-3

			double DampUsed = MIN2( 1e-2, t.Emis().damp() );
			double t1 = 3.02*pow(DampUsed*1e3,-0.064 );
			double u1 = sqrt(MAX2(tau,0.)*DampUsed )/SDIV(t1);
			wings = DampUsed/SDIV(t1)/sqrt( 0.78540 + POW2(u1) );

following botches result - some big

wolkje auto:cat serious.asr
agn_warm_absorber.out:  ChkMonitor botch>>IRON        9    -1.4061 =    -1.6370  -0.702   0.100  <<BIG BOTCH!!
agn_warm_absorber.out:  ChkMonitor botch>>IRON       10    -0.9514 =    -1.1490  -0.576   0.100  <<BIG BOTCH!!
agn_warm_absorber.out:  ChkMonitor botch>>IRON       11    -0.7946 =    -0.9535  -0.442   0.100  <<BIG BOTCH!!
agn_warm_absorber.out:  ChkMonitor botch>>IRON       12    -0.8528 =    -0.9557  -0.267   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       14    -1.2452 =    -1.1820   0.135   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       15    -1.4551 =    -1.3810   0.157   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       16    -1.5361 =    -1.4340   0.210   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       17    -0.7391 =    -0.5951   0.282   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       18    -0.9801 =    -0.8353   0.284   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       19    -1.5491 =    -1.4040   0.284   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       20    -2.6694 =    -2.5230   0.286   0.100
agn_warm_absorber.out:  ChkMonitor botch>>IRON       21    -4.1280 =    -3.9780   0.292   0.100
blr_level2.out:  ChkMonitor botch>>P  5 1117.98A     6.4717 =     6.4990   0.061   0.050
blr_level2.out:  ChkMonitor botch>>P  5 1128.01A     6.2315 =     6.2600   0.064   0.050
blr_level2.out:  ChkMonitor botch>>PHOS        5    -1.9042 =    -1.8690   0.078   0.050
blr_level2.out:  ChkMonitor botch>>PHOS        6    -2.9959 =    -2.9480   0.104   0.100
blr_n09_p20_Z20.out:  ChkMonitor botch>>totl 1035.00A     7.9767 =     8.0000   0.052   0.050
blr_n13_p18_Z20.out:  ChkMonitor botch>>totl 2326.00A     3.2033 =     3.2270   0.053   0.050
blr_nf84_45deg.out:  ChkMonitor botch>>HELI        3    -2.5043 =    -2.4820   0.050   0.050
blr_rnfa.out:  ChkMonitor botch>>c  3 977.000A     2.4194 =     2.5540   0.053   0.050
blr_rnfa.out:  ChkMonitor botch>>blnd 2798.00A    27.1653 =    25.5300  -0.064   0.050
blr_rnfa.out:  ChkMonitor botch>>Blnd 1888.00A     7.3824 =     7.7900   0.052   0.050
nlr_lex00.out:  ChkMonitor botch>>blnd 2424.00A     0.5035 =     0.5310   0.052   0.050
nlr_liner_grains.out:  ChkMonitor botch>>HELI        2    -0.7413 =    -0.7706  -0.070   0.050
nlr_paris.out:  ChkMonitor botch>>totl 1240.00A     0.2727 =     0.2979   0.085   0.050
nlr_paris.out:  ChkMonitor botch>>totl 1402.00A     0.5242 =     0.5596   0.063   0.050
nlr_paris.out:  ChkMonitor botch>>TOTL 1218.00A     0.2833 =     0.3132   0.095   0.050
nlr_paris.out:  ChkMonitor botch>>TOTL 1035.00A     0.1245 =     0.1452   0.143   0.100
nlr_paris.out:  ChkMonitor botch>>Blnd 1397.00A     0.2380 =     0.2527   0.058   0.050
pdr_leiden_f2.out:  ChkMonitor botch>>Si 2 34.8046m    -4.4268 =    -4.2150   0.386   0.050  <<BIG BOTCH!!
pdr_leiden_f4.out:  ChkMonitor botch>>Si 2 34.8046m    -4.3035 =    -4.2360   0.144   0.050
pdr_leiden_v2.out:  ChkMonitor botch>>Si 2 34.8046m    -3.5455 =    -3.5170   0.063   0.050
pn_fluc.out:  ChkMonitor botch>>blnd 2798.00A     0.9817 =     0.8152  -0.204   0.050  <<BIG BOTCH!!
pn_ots.out:  ChkMonitor botch>>blnd 2798.00A     1.2078 =     0.9973  -0.211   0.050  <<BIG BOTCH!!
pn_paris.out:  ChkMonitor botch>>blnd 2798.00A     1.2327 =     1.0590  -0.164   0.050  <<BIG BOTCH!!
pn_paris_cpre.out:  ChkMonitor botch>>blnd 2798.00A     1.3556 =     1.1280  -0.202   0.050  <<BIG BOTCH!!

comment:2 Changed 4 years ago by rjrw

comment:3 Changed 4 years ago by rjrw

Resolution: fixed
Status: newclosed

r8695 implements fixes to the Federman shielding function to bring it within -15% / +32% of the value obtained from exact integration across a wide range of tau, a (from -80% / +44% previously). As surmised in the original report, the min() does indeed maintain reasonable behaviour at small tau/large a.

With this change, the accuracy isn't substantially worse than the Rodgers & Williams fit used elsewhere in the literature (and now implemented as an option).

Note: See TracTickets for help on using tickets.