cloudy: /home66/gary/public_html/cloudy/c08_branch/source/cool

00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*CoolCalc compute calcium cooling */
00004 #include "cddefines.h"
00005 #include "taulines.h"
00006 #include "doppvel.h"
00007 #include "phycon.h"
00008 #include "ca.h"
00009 #include "dense.h"
00010 #include "thermal.h"
00011 #include "opacity.h"
00012 #include "rfield.h"
00013 #include "ligbar.h"
00014 #include "lines_service.h"
00015 #include "atoms.h"
00016 #include "cooling.h"
00017 
00018 void CoolCalc(void)
00019 {
00020         realnum p2;
00021         double  a21, 
00022           a31, 
00023           a41, 
00024           a42, 
00025           a51, 
00026           a52, 
00027           a53, 
00028           c21, 
00029           Ca2pop[5] ,
00030           cs, 
00031           cs2s2p, 
00032           cs2s3p        , 
00033           cs01, 
00034           cs02, 
00035           cs12, 
00036           cs14, 
00037           cs15, 
00038           d41, 
00039           d42, 
00040           d51, 
00041           d52, 
00042           d53, 
00043           hlgam, 
00044           op41, 
00045           op51, 
00046           opckh, 
00047           opcxyz, 
00048           PhotoRate2, 
00049           p3, 
00050           PhotoRate3, 
00051           PhotoRate4, 
00052           PhotoRate5, 
00053           r21, 
00054           r31, 
00055           r41, 
00056           r42, 
00057           r51, 
00058           r52, 
00059           r53;
00060         static double gCa2[5]={2.,4.,6.,2.,4.};
00061         static double exCa2[4]={13650.2,60.7,11480.6,222.9};
00062         static realnum opCax = 0.f;
00063         static realnum opCay = 0.f;
00064         static realnum opCaz = 0.f;
00065 
00066         DEBUG_ENTRY( "CoolCalc()" );
00067 
00068         /* Ca I 4228 */
00069         MakeCS(&TauLines[ipCaI4228]);
00070         atom_level2(&TauLines[ipCaI4228]);
00071 
00072         /* photoionization of evcited levels by Ly-alpha */
00073         hlgam = rfield.otslin[ Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont -1];
00074         PhotoRate5 = 1.7e-18*hlgam;
00075         PhotoRate4 = 8.4e-19*hlgam;
00076         PhotoRate3 = 7.0e-18*hlgam;
00077         PhotoRate2 = 4.8e-18*hlgam;
00078 
00079         /* spontaneous decays
00080          * frist two trans prob from 
00081          * >>refer      ca2     as      Zeippen, C.J. 1990, A&A, 229, 248 */
00082         a21 = 1.02*TauLines[ipT7324].Emis->Pesc;
00083         a31 = 1.05*TauLines[ipT7291].Emis->Pesc;
00084         a41 = 1.4e8*TauLines[ipT3969].Emis->Pesc;
00085         a51 = 1.4e8*TauLines[ipT3934].Emis->Pesc;
00086         a42 = 7.9e6*TauLines[ipT8662].Emis->Pesc;
00087         a52 = 8.2e5*TauLines[ipT8498].Emis->Pesc;
00088         a53 = 7.48e6*TauLines[ipT8542].Emis->Pesc;
00089 
00090         /* destruction of IR triplet by continuous opacities */
00091         opcxyz = opac.opacity_abs[ TauLines[ipT7324].ipCont -1];
00092 
00093         /* opcxyz = opac(icaxyz) */
00094         if( opcxyz > 0. )
00095         {
00096                 d52 = 5.6*opcxyz/(opcxyz + opCax)*(1. - TauLines[ipT8498].Emis->Pesc);
00097                 d53 = 5.6*opcxyz/(opcxyz + opCay)*(1. - TauLines[ipT8542].Emis->Pesc);
00098                 d42 = 5.6*opcxyz/(opcxyz + opCaz)*(1. - TauLines[ipT8662].Emis->Pesc);
00099         }
00100         else
00101         {
00102                 d52 = 0.;
00103                 d53 = 0.;
00104                 d42 = 0.;
00105         }
00106 
00107         /* near UV dest of KH by background continuum */
00108         opckh = opac.opacity_abs[ TauLines[ipT3969].ipCont -1];
00109 
00110         /* opckh = opac(icakh) */
00111         if( opckh > 0. )
00112         {
00113                 op51 = dense.xIonDense[ipCALCIUM][1]*3.89e-7/DoppVel.doppler[19];
00114                 d51 = 5.6*opckh/(opckh + op51);
00115                 op41 = dense.xIonDense[ipCALCIUM][1]*1.96e-7/DoppVel.doppler[19];
00116                 d41 = 5.6*opckh/(opckh + op41);
00117         }
00118         else
00119         {
00120                 op51 = 0.;
00121                 d51 = 0.;
00122                 op41 = 0.;
00123                 d41 = 0.;
00124         }
00125         /* net rates */
00126         r21 = PhotoRate2 + a21;
00127         r31 = PhotoRate3 + a31;
00128         r41 = a41 + PhotoRate4 + d41;
00129         r51 = a51 + PhotoRate5 + d51;
00130         r42 = a42 + d42;
00131         r52 = a52 + d52;
00132         r53 = a53 + d53;
00133         cs14 = 0.923*phycon.te10*phycon.te10;
00134         cs15 = cs14*2.;
00135         TauLines[ipT3969].Coll.cs = (realnum)cs14;
00136         TauLines[ipT3934].Coll.cs = (realnum)cs15;
00137 
00138         /* following used to correct rec contribution
00139          * fcakh = a51 / ( a51 + eden*1.5e-5 / sqrte ) 
00140          * cs 1-2 from 
00141          * >>refer      ca2     cs      Saraph, H.E. 1970, J.Phys. B, 3, 952
00142          * other 
00143          * >>refer      ca2     cs      Chidichimo, M.C. 1981, J.Phys. B, 14, 4149 */
00144         atom_pop5(gCa2,exCa2,5.8,8.6,cs14,cs15,20.6,22.9,9.8,3.4,44.4,1.0,
00145           r21,r31,r41,r51,0.,r42,r52,0.,r53,0.,Ca2pop,dense.xIonDense[ipCALCIUM][1]);
00146 
00147         /* CDSQTE = 8.629E-6*EDEN/SQRTE */
00148         c21 = 5.8/4.*dense.cdsqte;
00149 
00150         /* remember largest ratio of Ly-al removal to total */
00151         if( dense.xIonDense[ipCALCIUM][1] > 0. )
00152         {
00153                 ca.Ca2RmLya = MAX2(ca.Ca2RmLya,(realnum)(PhotoRate2/(PhotoRate2+a21+c21)));
00154         }
00155         ca.Cak = (realnum)(Ca2pop[4]*a51*5.06e-12);
00156         ca.Cah = (realnum)(Ca2pop[3]*a41*5.01e-12);
00157         ca.Cax = (realnum)(Ca2pop[4]*a52*2.34e-12);
00158         ca.Cay = (realnum)(Ca2pop[4]*a53*2.33e-12);
00159         ca.Caz = (realnum)(Ca2pop[3]*a42*2.30e-12);
00160         ca.Caf1 = (realnum)(Ca2pop[2]*a31*2.73e-12);
00161         ca.Caf2 = (realnum)(Ca2pop[1]*a21*2.72e-12);
00162         ca.popca2ex = (realnum)(Ca2pop[1] + Ca2pop[2] + Ca2pop[3] + Ca2pop[4]);
00163 
00164         /* this is the total cooling due to the model atom */
00165         ca.Cair = ca.Cax + ca.Cay + ca.Caz;
00166         ca.c7306 = ca.Caf1 + ca.Caf2;
00167         ca.Cakh = ca.Cak + ca.Cah;
00168 
00169         /* add the CaII lines to the cooling stack */
00170         /* >>chng 04 feb 28, next three, deriv wrt temp had not been present 
00171          * extra factor of three in second term of deriv is due to extra T-1
00172          * dependence of cooling on recombination coefficient.  Almost always
00173          * CaII cooling is most important when Ca+ is a trace stage of ionization
00174          * so changes in temp change cooling due to change in Ca+ abundance */
00175         CoolAdd("Ca 2",7306,ca.c7306);
00176         thermal.dCooldT -= ca.c7306*4.*thermal.halfte;
00177 
00178         CoolAdd("Ca 2",8400,ca.Cair);
00179         /* >>chng 04 feb 28, this is emperical fit to cooling deriv - does down with
00180          * increasing temp */
00181         thermal.dCooldT -= ca.Cair*4.*thermal.halfte;
00182 
00183         /* >>chng 04 feb 28 wl had been 3800, chng to 3950, closer to mean of K and R */
00184         CoolAdd("Ca 2",3950,ca.Cakh);
00185         /* >>chng 04 feb 28, this is emperical fit to cooling deriv - does down with
00186          * increasing temp */
00187         thermal.dCooldT -= ca.Cakh*4.*thermal.halfte;
00188 
00189         /*fprintf(ioQQQ,"DEBUG ca2\t%.2f\t%.5e\t%.4e\t%.4e\n",
00190                 fnzone, phycon.te,ca.Cakh,dense.xIonDense[ipCALCIUM][1]);*/
00191 
00192         /* level populations that will be used for excited state photoionization */
00193         ca.dstCala = (realnum)(Ca2pop[4]*PhotoRate5 + Ca2pop[3]*PhotoRate4);
00194         ca.dCakh = (realnum)(ca.dstCala*5.03e-12);
00195         ca.dCaf12 = (realnum)((Ca2pop[2]*PhotoRate3 + Ca2pop[1]*PhotoRate2)*2.31e-12);
00196         opCax = (realnum)(Ca2pop[1]*1.13e-8/DoppVel.doppler[19]);
00197         opCay = (realnum)(Ca2pop[2]*6.87e-8/DoppVel.doppler[19]);
00198         opCaz = (realnum)(Ca2pop[1]*5.74e-8/DoppVel.doppler[19]);
00199 
00200         /* total rate Lalpha destroys CaII,
00201          * this is only used in ioncali to increase ionization rate by
00202          * adding it to the ct vector */
00203         if( dense.xIonDense[ipCALCIUM][1] > 0. )
00204         {
00205                 ca.dstCala = (realnum)(
00206                         (ca.dstCala + ca.dCaf12/2.31e-12)/dense.xIonDense[ipCALCIUM][1]);
00207                 /* take average of old and new */
00208                 ca.dstCala = (realnum)((ca.dstCala + ca.oldcla)/2.);
00209                 ca.oldcla = ca.dstCala;
00210                 {
00211                         /*@-redef@*/
00212                         enum {DEBUG_LOC=false};
00213                         /*@+redef@*/
00214                         if( DEBUG_LOC )
00215                         {
00216                                 fprintf(ioQQQ," hlgam is %e\n", hlgam);
00217                         }
00218                 }
00219         }
00220         else
00221         {
00222                 ca.dstCala = 0.;
00223         }
00224         ca.Ca3d = (realnum)(Ca2pop[1] + Ca2pop[2]);
00225         ca.Ca4p = (realnum)(Ca2pop[3] + Ca2pop[4]);
00226 
00227         /* incl stimulated emission for Calcium II 5-level atom */
00228         TauLines[ipT3934].Emis->PopOpc = (Ca2pop[0] - Ca2pop[4]/2.);
00229         TauLines[ipT3934].Hi->Pop = Ca2pop[4];
00230         TauLines[ipT3934].Lo->Pop = Ca2pop[0];
00231         TauLines[ipT3969].Emis->PopOpc = (Ca2pop[0] - Ca2pop[3]);
00232         TauLines[ipT3969].Hi->Pop = Ca2pop[3];
00233         TauLines[ipT3969].Lo->Pop = Ca2pop[0];
00234 
00235         opac.tpcah[0] = TauLines[ipT3969].Emis->TauIn;
00236         TauLines[ipT8498].Emis->PopOpc = (Ca2pop[1] - Ca2pop[4]);
00237         TauLines[ipT8498].Hi->Pop = Ca2pop[4];
00238         TauLines[ipT8498].Lo->Pop = Ca2pop[1];
00239         TauLines[ipT8542].Emis->PopOpc = (Ca2pop[2] - Ca2pop[4]*1.5);
00240         TauLines[ipT8542].Hi->Pop = Ca2pop[4];
00241         TauLines[ipT8542].Lo->Pop = Ca2pop[2];
00242         TauLines[ipT8662].Emis->PopOpc = (Ca2pop[1] - Ca2pop[3]*2.);
00243         TauLines[ipT8662].Hi->Pop = Ca2pop[3];
00244         TauLines[ipT8662].Lo->Pop = Ca2pop[1];
00245         TauLines[ipT7291].Emis->PopOpc = dense.xIonDense[ipCALCIUM][1];
00246         TauLines[ipT7291].Hi->Pop = 0.;
00247         TauLines[ipT7291].Lo->Pop = dense.xIonDense[ipCALCIUM][1];
00248         TauLines[ipT7324].Emis->PopOpc = dense.xIonDense[ipCALCIUM][1];
00249         TauLines[ipT7324].Hi->Pop = 0.;
00250         TauLines[ipT7324].Lo->Pop = dense.xIonDense[ipCALCIUM][1];
00251 
00252         /* Ca IV 3.2 micron; data from
00253          * >>refer      ca4     as      Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
00254          * >>refercon   ed by D.R. Flower, (D. Reidel: Holland), 143
00255          * Y(ik) from 
00256          * >>refer      ca4     cs      Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209 */
00257         if( phycon.te <= 1e5 )
00258         {
00259                 cs = MAX2(1.0,8.854e-3*phycon.sqrte);
00260         }
00261         else if( phycon.te < 2.512e5 )
00262         {
00263                 cs = 2.8;
00264         }
00265         else
00266         {
00267                 cs = 641.1/(phycon.te30*phycon.te10*phycon.te02*phycon.te02/
00268                   phycon.te003);
00269         }
00270         PutCS(cs,&TauLines[ipTCa3]);
00271         atom_level2(&TauLines[ipTCa3]);
00272 
00273         /* [Ca V] IR 4.16, 11.47 micron; A from
00274          * >>refer      ca5     as      Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
00275          * >>refercon   ed by D.R. Flower, (D. Reidel: Holland), 143
00276          * cs from 
00277          * >>refer      ca5     cs      Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
00278          * >>chng 96 jul 16, big changes in cs */
00279         cs = MIN2(3.3,0.392*phycon.te20/phycon.te005/phycon.te003);
00280         cs = MAX2(2.2,cs);
00281         PutCS(cs,&TauLines[ipTCa4]);
00282 
00283         /* >>chng 96 aug 02, following had error in te dep in ver 90.01 */
00284         cs = MIN2(0.93,0.162*phycon.te10*phycon.te05*phycon.te003*
00285           phycon.te001);
00286         cs = MAX2(0.67,cs);
00287         PutCS(cs,&TauLines[ipTCa12]);
00288 
00289         cs = MIN2(0.97,0.0894*phycon.te20*phycon.te01*phycon.te005);
00290         cs = MAX2(0.60,cs);
00291         PutCS(cs,&TauDummy);
00292 
00293         atom_level3(&TauLines[ipTCa4],&TauLines[ipTCa12],&TauDummy);
00294 
00295         /* Ca V lines from 1d, 1s; A from 
00296          * >>refer      ca5     as      Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
00297          * >>refercon   ed by D.R. Flower, (D. Reidel: Holland), 143
00298          * cs from 
00299          * >>refer      ca5     cs       Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
00300          * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
00301         cs01 = MIN2(4.1,0.533*phycon.te20/phycon.te01);
00302         cs01 = MAX2(2.8,cs01);
00303         cs02 = MIN2(0.87,5.22e-03*phycon.sqrte);
00304         p3 = atom_pop3(9.,5.,1.,cs01,cs02,1.35,2.326,23.2,3.73,2.57e4,3.60e4,
00305           &p2,dense.xIonDense[ipCALCIUM][4],0.,0.,0.);
00306 
00307         ca.c3997 = p3*3.73*4.98e-12;
00308         ca.c2414 = p3*23.1*8.245e-12;
00309         ca.Ca6087 = p2*0.426*3.268e-12;
00310         ca.c5311 = p2*1.90*3.747e-12;
00311 
00312         CoolAdd("Ca 5",3997,ca.c3997);
00313         CoolAdd("Ca 5",2414,ca.c2414);
00314         CoolAdd("Ca 5",6087,ca.Ca6087);
00315         CoolAdd("Ca 5",5311,ca.c5311);
00316 
00317         /* Ca VII lines from 1d, 1s
00318          * all cs from 
00319          * >>refer      ca72    cs      Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
00320          * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
00321         cs01 = MIN2(4.4,22.25/(phycon.te20/phycon.te02/phycon.te02));
00322         cs01 = MAX2(3.5,cs01);
00323 
00324         cs12 = MIN2(1.20,0.303*phycon.te30*phycon.te03);
00325         cs12 = MAX2(0.62,cs12);
00326 
00327         cs02 = MIN2(0.959,7.889/(phycon.te20*phycon.te05/phycon.te01));
00328         cs02 = MAX2(0.50,cs02);
00329 
00330         p3 = atom_pop3(9.,5.,1.,cs01,cs02,cs12,3.124,30.4,6.81,2.91e4,3.90e4,
00331           &p2,dense.xIonDense[ipCALCIUM][6],0.,0.,0.);
00332 
00333         ca.Ca3688 = p3*6.81*5.40e-12;
00334         ca.Ca2112 = p3*30.4*9.42e-12;
00335         ca.Ca5620 = p2*2.15*3.548e-12;
00336         ca.Ca4941 = p2*0.974*4.037e-12;
00337         CoolAdd("Ca 7",3688,ca.Ca3688);
00338         CoolAdd("Ca 7",2112,ca.Ca2112);
00339         CoolAdd("Ca 7",5620,ca.Ca5620);
00340         CoolAdd("Ca 7",4941,ca.Ca4941);
00341 
00342         /* all cs from 
00343          * >>refer      ca7     cs      Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
00344          * [Ca VII] 4.09, 6.15 mic 3P lines */
00345         cs = MIN2(5.354,0.406*phycon.te20*phycon.te03*phycon.te01);
00346         cs = MAX2(3.702,cs);
00347         PutCS(cs,&TauLines[ipCa0741]);
00348 
00349         cs = MIN2(1.59,0.183*phycon.te20);
00350         cs = MAX2(1.153,cs);
00351         PutCS(cs,&TauLines[ipCa0761]);
00352 
00353         cs = MIN2(1.497,0.0917*phycon.te20*phycon.te05* phycon.te01);
00354         cs = MAX2(1.005,cs);
00355         PutCS(cs,&TauDummy);
00356 
00357         /* atom_level3(  t10,t21,t20) */
00358         atom_level3(&TauLines[ipCa0761],&TauLines[ipCa0741],&TauDummy);
00359 
00360         /* [Ca VIII]  2.32 microns, cs 
00361          * >>refer      ca8     cs      Saraph, H.E., & Storey, P.J. A&AS, 115, 151 */
00362         cs = MIN2(6.75,22.04/(phycon.te10*phycon.te02*phycon.te005));
00363         PutCS(cs,&TauLines[ipCa08232]);
00364         atom_level2(&TauLines[ipCa08232]);
00365 
00366         /* [Ca 12] 3328.78A, cs from 
00367          * >>refer      ca12    cs      Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
00368         cs = MIN2(0.172,0.0118*phycon.te20*phycon.te01);
00369         cs = MAX2(0.10,cs);
00370         PutCS(cs,&TauLines[ipCa12333]);
00371         atom_level2(&TauLines[ipCa12333]);
00372 
00373         /* Li seq Ca 18 2s2p 2s3p, 2s2p as two separate lines
00374          * >>refer      ca18    cs      Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
00375         ligbar(20,&TauLines[ipTCa302],&TauLines[ipTCa19],&cs2s2p,&cs2s3p);
00376 
00377         PutCS(cs2s2p,&TauLines[ipTCa302]);
00378         atom_level2(&TauLines[ipTCa302]);
00379 
00380         /* funny factor (should have been 0.5) due to energy change */
00381         PutCS(cs2s2p*0.439,&TauLines[ipTCa345]);
00382         atom_level2(&TauLines[ipTCa345]);
00383 
00384         PutCS(cs2s3p,&TauLines[ipTCa19]);
00385         atom_level2(&TauLines[ipTCa19]);
00386         return;
00387 }
/home66/gary/public_html/cloudy/c08_branch/source/cool_calc.cpp
