Changeset 1917
- Timestamp:
- 04/05/08 17:25:20 (9 months ago)
- Location:
- trunk/source
- Files:
-
- 2 modified
-
mole_h2_coll.cpp (modified) (9 diffs)
-
prt_met.cpp (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/mole_h2_coll.cpp
r1815 r1917 21 21 static bool **lgDefn_H2He_coll; 22 22 23 STATIC realnum H2_CollidRateEvalOne( long iVibHi, long iRotHi,long iVibLo, 24 long iRotLo, long ipHi , long ipLo , long nColl ) 23 /* compute rate coefficient for a single quenching collision */ 24 STATIC realnum H2_CollidRateEvalOne( 25 /*returns collision rate coefficient, cm-3 s-1 for quenching collision 26 * from this upper state */ 27 long iVibHi, long iRotHi,long iVibLo, 28 /* to this lower state */ 29 long iRotLo, long ipHi , long ipLo , 30 /* colliders are H, He, H2(ortho), H2(para), and H+ */ 31 long nColl ) 25 32 { 26 33 double fitted; 27 34 realnum rate; 28 double t = phycon.te/1000. + 1.;29 double t 2 = POW2(t);35 double t3Plus1 = phycon.te/1000. + 1.; 36 double t3Plus1Squared = POW2(t3Plus1); 30 37 /* these are fits to the existing collision data 31 38 * used to create g-bar rates */ … … 68 75 } 69 76 70 /* >>chng 07 apr 2, GS 71 * this is option to use g-bar g bar guess of rate coefficient for 72 * collision with He, this does not change ortho & para 73 * code had always applied g-bar, which was incorrect 74 * turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 77 /* use g-bar g bar guess of rate coefficient for 78 * collision with He, this does not change ortho & para 79 * turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 75 80 else if( mole.lgColl_gbar && 76 81 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) ) … … 139 144 else if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0]!= 0 ) 140 145 { 141 /* evaluate collision rates for those with real collision data */ 146 /* these are the fits from 147 *>>refer H2 coll Le Bourlot, J., Pineau des Forets, 148 *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802 149 * evaluate collision rates for those with real collision data */ 142 150 143 151 double r = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] + 144 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t +145 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t 2;152 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t3Plus1 + 153 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t3Plus1Squared; 146 154 147 155 rate = (realnum)pow(10.,r)*mole.lgColl_deexec_Calc; … … 201 209 } 202 210 203 /*H2_CollidRateEvalAll - set H2 collision rate s */211 /*H2_CollidRateEvalAll - set H2 collision rate coefficients */ 204 212 void H2_CollidRateEvalAll( void ) 205 213 { … … 302 310 303 311 /* >>chng 02 nov 13, the Sun and Dalgarno rates diverge to + inf below this temp */ 304 /* >>chng 05 feb 09, do not return zero when t< 100 - instead, don't let T fall below 100 */312 /* >>chng 05 feb 09, do not return zero when T < 100 - instead, don't let T fall below 100 */ 305 313 double excit1; 306 314 double te_used = MAX2( 100. , phycon.te ); … … 657 665 { 658 666 659 double k, t , b2, c2, f2, t2;667 double k, t3Plus1, b2, c2, f2, t3; 660 668 661 669 DEBUG_ENTRY( "H2_He_coll()" ); … … 695 703 temp = MIN2( 1e4 , temp ); 696 704 697 t 2= temp/1e3;698 /* t = T*10^-3+ 1*/699 t = t2+1;700 b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t 2+1);701 c2 = H2_He_coll_fit_par[init][final][2]/(t *t);705 t3 = temp/1e3; 706 /* t3Plus1 = T/1000 + 1*/ 707 t3Plus1 = t3+1; 708 b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t3+1); 709 c2 = H2_He_coll_fit_par[init][final][2]/(t3Plus1*t3Plus1); 702 710 { 703 711 enum {DEBUG_LOC=false}; … … 707 715 init,final, 708 716 H2_He_coll_fit_par[init][final][5], 709 H2_He_coll_fit_par[init][final][6]*t 2+1,717 H2_He_coll_fit_par[init][final][6]*t3+1, 710 718 H2_He_coll_fit_par[init][final][7] 711 /*pow(H2_He_coll_fit_par[init][final][6]*t 2+1,H2_He_coll_fit_par[init][final][7])*/719 /*pow(H2_He_coll_fit_par[init][final][6]*t3+1,H2_He_coll_fit_par[init][final][7])*/ 712 720 ); 713 721 } 714 722 } 715 723 /* this is log of f2 - see whether it is within bounds */ 716 sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t 2+1. );724 sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t3+1. ); 717 725 /* this protects against overflow */ 718 726 if( fabs(sum1)< 38. ) … … 721 729 * bug in the pgcc compiler that caused h2_pdr_leiden_f1 to throw an fpe 722 730 * the changed code is equivalent */ 723 /*f2 = H2_He_coll_fit_par[init][final][5]/pow(H2_He_coll_fit_par[init][final][6]*t 2+1.,H2_He_coll_fit_par[init][final][7]);*/731 /*f2 = H2_He_coll_fit_par[init][final][5]/pow(H2_He_coll_fit_par[init][final][6]*t3+1.,H2_He_coll_fit_par[init][final][7]);*/ 724 732 f2 = H2_He_coll_fit_par[init][final][5]/pow(10. , sum1 ); 725 733 } -
trunk/source/prt_met.cpp
r1914 r1917 25 25 if( prt.lgPrtTau ) 26 26 { 27 fprintf( ioQQQ, " Line Optical Depths\n"); 27 28 28 /* 'In' means to zero counter, following elements mean nothing*/29 prme("I n",&TauLines[0]);29 /* "IN" - initialize */ 30 prme("IN",&TauLines[0]); 30 31 31 /* do hydrogenic species */ 32 /* >>chng 02 may 16, bring two loops together */ 32 /* iso sequences */ 33 33 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 34 34 { … … 54 54 prme(" c",&TauLines[i]); 55 55 } 56 57 /*Atomic & Molecular Lines*/58 printf("\n Atomic and Molecular Lines \n");59 for( i=0; i < linesAdded2; i++)60 {61 prme(" am",atmolEmis[i].tran);62 }63 printf(" \n");64 56 65 57 for( i=0; i < nWindLine; i++ ) … … 95 87 } 96 88 97 /* this means print it */ 98 prme(" ",&TauLines[0]); 89 /* data base lines */ 90 for( i=0; i < linesAdded2; i++) 91 { 92 prme("DB",atmolEmis[i].tran); 93 } 94 95 fprintf( ioQQQ, "\n"); 99 96 } 100 97 return; 101 98 } 102 99 103 /*prme print heavy element line optical depths at end of calculation */ 104 #define NLIM 6 105 106 /****************************************************************************** 107 * 108 * prme - the helper routine for above 109 * 110 ******************************************************************************/ 111 100 /* prme - print optical depth */ 112 101 void prme( 102 /* flag saying what to do 103 * "IN" initialize 104 * " c" add to list of old style lines 105 * "DB" add to list of database lines 106 */ 113 107 const char *chDoIt, 114 108 transition * t) 115 109 { 116 static char chLsav[NLIM][11]; 117 //char chAtMolWL[20],chAtMol[35]; 118 long int i; 119 static long int n; 120 /*lint -e785 too few init for aggregate */ 121 static realnum tsav[NLIM] = {0.}; 122 /*lint +e785 too few init for aggregate */ 110 char chAtMolWL[20],chAtMol[35]; 111 static long int n ; 123 112 124 113 DEBUG_ENTRY( "prme()" ); 125 114 126 /* this routine is called by prtmet, and is used to print out127 * line optical depths after model complete - it IS NOT a128 * timing pace setter */129 130 /* label for line */131 /* strcpy( chLabel, elementnames.chElementSym[ t->nelem -1] ); */132 133 115 if( t->ipCont <= 0 ) 134 116 { 117 /* line is not transferred */ 135 118 return; 136 119 } 137 120 138 if( strcmp(chDoIt,"In") == 0 ) 121 /* "In" is to initialize for new printout */ 122 if( strncmp(chDoIt,"IN",2) == 0 ) 139 123 { 140 124 n = 0; 125 } 141 126 127 else if( strncmp(chDoIt,"DB",2) == 0) 128 { 129 /* database lines, - cannot now address species labels for atoms 130 * and molecules in one simple way so separate from most lines 131 * print optical depth if greater than lower limit, or significantly negative */ 132 if( t->Emis->TauIn > prt.PrtTauFnt || t->Emis->TauIn < -1e-5 ) 133 { 134 sprt_wl(chAtMolWL,t->WLAng); 135 strcpy(chAtMol,t->Hi->chLabel); 136 strcat(chAtMol," "); 137 strcat(chAtMol,chAtMolWL); 138 fprintf( ioQQQ, " %10.15s",chAtMol); 139 fprintf( ioQQQ,PrintEfmt("%9.2e", t->Emis->TauIn)); 140 fprintf( ioQQQ, " "); 141 // throw CR after printing 6 numbers 142 ++n; 143 if(n == 6) 144 { 145 n = 0; 146 fprintf( ioQQQ, " \n"); 147 } 148 } 142 149 } 143 else if( strcmp(chDoIt," ") == 0 ) 150 151 else if( strncmp(chDoIt," c",2) == 0) 144 152 { 145 /* NB number in following statement must agree with nlim */ 146 fprintf( ioQQQ, " " ); 147 for( i=0; i < n; i++ ) 153 /* print optical depth if greater than lower limit, or significantly negative */ 154 if( t->Emis->TauIn > prt.PrtTauFnt || t->Emis->TauIn < -1e-5 ) 148 155 { 149 fprintf( ioQQQ, "%10.10s",chLsav[i]); 150 fprintf( ioQQQ,PrintEfmt("%9.2e", tsav[i] )); 151 fprintf( ioQQQ, " "); 156 /* PrtTauFnt is threshold for printing it */ 157 fprintf( ioQQQ, " %10.10s",chLineLbl(t)); 158 fprintf( ioQQQ, PrintEfmt("%9.2e", t->Emis->TauIn )); 159 160 // throw CR after printing 6 numbers 161 ++n; 162 if(n == 6) 163 { 164 n = 0; 165 fprintf( ioQQQ, " \n"); 166 } 152 167 } 153 fprintf( ioQQQ, " \n" );154 155 168 } 156 /* print optical depth if greater than lower limit, or significantly negative */ 157 else if( t->Emis->TauIn > prt.PrtTauFnt || t->Emis->TauIn < -1e-5 ) 158 { 159 /* t->t(ipLnTauIn) is the optical depth 160 * PrtTauFnt is threshold for printing it */ 161 if( n >= NLIM ) 162 { 163 /* NB number in following statement must agree with nlim */ 164 fprintf( ioQQQ, " " ); 165 for( i=0; i < NLIM; i++ ) 166 { 167 fprintf( ioQQQ, "%10.10s",chLsav[i]); 168 fprintf( ioQQQ,PrintEfmt("%9.2e", tsav[i] )); 169 fprintf( ioQQQ, " "); 170 } 171 fprintf( ioQQQ, " \n" ); 172 n = 0; 173 } 174 175 strcpy( chLsav[n], chLineLbl(t) ); 176 tsav[n] = t->Emis->TauIn; 177 ++n; 178 } 169 else 170 TotalInsanity(); 179 171 return; 180 172 }
