cloudy trunk
|
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 /*PunchLineData punches selected line data for all lines transferred in code */ 00004 /*Punch1LineData punch data for one line */ 00005 #include "cddefines.h" 00006 #include "taulines.h" 00007 #include "hmi.h" 00008 #include "iso.h" 00009 #include "phycon.h" 00010 #include "physconst.h" 00011 #include "elementnames.h" 00012 #include "hydrogenic.h" 00013 #include "lines_service.h" 00014 #include "dense.h" 00015 #include "atomfeii.h" 00016 #include "lines.h" 00017 #include "atmdat.h" 00018 #include "prt.h" 00019 #include "mole_co_atom.h" 00020 #include "h2.h" 00021 #include "thermal.h" 00022 #include "cooling.h" 00023 #include "punch.h" 00024 00025 NORETURN void PunchLineData(FILE * ioPUN) 00026 { 00027 00028 long int i, 00029 j, 00030 limit , 00031 nelem , 00032 ipHi , 00033 ipLo; 00034 00035 const long nskip=2; /* number of emission lines per line of output */ 00036 double tot; 00037 bool lgElemOff=false; 00038 realnum a , b; /* dummy vars to pass to rotate cooling routine */ 00039 00040 DEBUG_ENTRY( "PunchLineData()" ); 00041 00042 /* routine punches out (on unit ioPUN) line data 00043 * for all recombination lines, and all transitions that are transferred */ 00044 00045 /* say what is happening so we know why we stopped */ 00046 fprintf( ioQQQ, " punching line data, then stopping\n" ); 00047 00048 /* first check that all lines are turned on */ 00049 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00050 { 00051 if( !dense.lgElmtOn[nelem] ) 00052 { 00053 fprintf(ioQQQ," WARNING - I am punching line data but element %s is turned off.\n", 00054 elementnames.chElementName[nelem]); 00055 lgElemOff = true; 00056 } 00057 } 00058 if( lgElemOff ) 00059 { 00060 fprintf(ioQQQ,"Some elements are turned off and punch line data requested.\n"); 00061 fprintf(ioQQQ,"Code is now designed to do punch line data only with all elements on.\n"); 00062 fprintf(ioQQQ,"Please try again with all elements on.\n"); 00063 cdEXIT(EXIT_FAILURE); 00064 } 00065 00066 /* evaluate rec coefficient at constant temperature if this is set, else 00067 * use 10,000K */ 00068 double TeNew; 00069 if( thermal.lgTemperatureConstant ) 00070 { 00071 TeNew = thermal.ConstTemp; 00072 } 00073 else 00074 { 00075 TeNew = 1e4; 00076 } 00077 TempChange(TeNew , false); 00078 00079 /* this is set of Dima's recombination lines */ 00080 t_ADfA::Inst().rec_lines(phycon.te,LineSave.RecCoefCNO); 00081 fprintf( ioPUN, "\n Recombination lines of C, N, O\n" ); 00082 fprintf( ioPUN, " Ion WL(A) Coef Ion WL(A) Coef\n" ); 00083 for( i=0; i<471; i+=nskip) 00084 { 00085 /* nskip is set to 2 above */ 00086 limit = MIN2(471,i+nskip); 00087 fprintf( ioPUN, " " ); 00088 for( j=i; j < limit; j++ ) 00089 { 00090 fprintf( ioPUN, "%2.2s%2ld%6ld%8.3f ", 00091 elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][j])-1], 00092 (long)(LineSave.RecCoefCNO[0][j]-LineSave.RecCoefCNO[1][j]+1.01), 00093 (long)(LineSave.RecCoefCNO[2][j]+0.5), 00094 log10(SDIV(LineSave.RecCoefCNO[3][j]) ) ); 00095 } 00096 fprintf( ioPUN, " \n" ); 00097 } 00098 fprintf( ioPUN, "\n\n" ); 00099 00100 dense.eden = 1.; 00101 dense.gas_phase[ipHYDROGEN] = 1.; 00102 dense.EdenHCorr = 1.; 00103 00104 /* want very small neutral fractions so get mostly e- cs */ 00105 dense.xIonDense[ipHYDROGEN][1] = 1.e-5f; 00106 hmi.Hmolec[ipMH2g] = 0.; 00107 dense.xIonDense[ipHYDROGEN][1] = 1.; 00108 for( i=1; i <= nLevel1; i++ ) 00109 { 00110 TauLines[i].Lo->Pop = 1.; 00111 } 00112 00113 for( i=0; i < nWindLine; i++ ) 00114 { 00115 TauLine2[i].Lo->Pop = 1.; 00116 } 00117 00118 for( i=0; i < nUTA; i++ ) 00119 { 00120 UTALines[i].Lo->Pop = 1.; 00121 } 00122 00123 for( i=0; i < LIMELM; i++ ) 00124 { 00125 for( j=0; j < LIMELM+1; j++ ) 00126 { 00127 dense.xIonDense[i][j] = 1.; 00128 } 00129 } 00130 00131 /* evaluate cooling, this forces evaluation of collision strengths */ 00132 CoolEvaluate(&tot); 00133 00134 fprintf( ioPUN, " Level 1 transferred lines\n" ); 00135 00136 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00137 00138 for( i=1; i <= nLevel1; i++ ) 00139 { 00140 /* chLineLbl generates a 1 char string from the line transfer array info - 00141 * "Ne 2 128" the string is null terminated, 00142 * in following printout the first 4 char is used first, followed by 00143 * an integer, followed by the rest of the array*/ 00144 Punch1LineData( &TauLines[i] , ioPUN , true); 00145 } 00146 00147 fprintf( ioPUN, "\n\n\n" ); 00148 fprintf( ioPUN, " end level 1, start level 2\n" ); 00149 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00150 for( i=0; i < nWindLine; i++ ) 00151 { 00152 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00153 { 00154 Punch1LineData( &TauLine2[i] , ioPUN , true); 00155 } 00156 } 00157 00158 fprintf( ioPUN, "\n\n\n" ); 00159 fprintf( ioPUN, " end level 2, start inner shell\n" ); 00160 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00161 00162 for( i=0; i < nUTA; i++ ) 00163 { 00164 if( UTALines[i].Emis->Aul > 0. ) 00165 Punch1LineData( &UTALines[i] , ioPUN , true); 00166 } 00167 00168 fprintf( ioPUN, "\n\n\n" ); 00169 fprintf( ioPUN, " end inner shell, start h-like iso seq\n" ); 00170 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00171 00172 /* h-like iso sequence */ 00173 /* the hydrogen like iso-sequence */ 00174 for( nelem=0; nelem < LIMELM; nelem++ ) 00175 { 00176 iso_collide( ipH_LIKE, nelem ); 00177 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00178 { 00179 /* arrays are dim'd iso.numLevels_max[ipH_LIKE][nelem]+1 */ 00180 /* keep this limit to iso.numLevels_max, instead of _local. */ 00181 for( ipLo=ipH1s; ipLo < iso.numLevels_max[ipH_LIKE][nelem]-1; ipLo++ ) 00182 { 00183 for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ ) 00184 { 00185 Punch1LineData( &Transitions[ipH_LIKE][nelem][ipHi][ipLo] , ioPUN , false ); 00186 } 00187 } 00188 } 00189 } 00190 00191 fprintf( ioPUN, "\n\n\n" ); 00192 fprintf( ioPUN, " end h-like iso seq, start he-like iso seq\n" ); 00193 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00194 for( nelem=1; nelem < LIMELM; nelem++ ) 00195 { 00196 if( nelem < 2 || dense.lgElmtOn[nelem] ) 00197 { 00198 /* arrays are dim'd iso.numLevels_max[ipH_LIKE][nelem]+1 */ 00199 for( ipLo=ipHe1s1S; ipLo < iso.numLevels_max[ipHE_LIKE][nelem]-1; ipLo++ ) 00200 { 00201 for( ipHi=ipLo+1; ipHi < iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ ) 00202 { 00203 Punch1LineData( &Transitions[ipHE_LIKE][nelem][ipHi][ipLo] , ioPUN , false ); 00204 } 00205 } 00206 } 00207 } 00208 00209 fprintf( ioPUN, "\n\n\n" ); 00210 fprintf( ioPUN, " end he-like iso seq, start hyperfine structure lines\n" ); 00211 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00212 /* fine structure lines */ 00213 for( i=0; i < nHFLines; i++ ) 00214 { 00215 Punch1LineData( &HFLines[i] , ioPUN , true); 00216 } 00217 00218 fprintf( ioPUN, "\n\n\n" ); 00219 fprintf( ioPUN, " end hyperfine, start database lines\n" ); 00220 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00221 /* Databases: Atoms & Molecules*/ 00222 for( i=0; i < linesAdded2; i++ ) 00223 { 00224 Punch1LineData( atmolEmis[i].tran , ioPUN , true); 00225 } 00226 00227 for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ ) 00228 { 00229 for( long nelem = ipISO; nelem < LIMELM; nelem++ ) 00230 { 00231 /* must always do helium even if turned off */ 00232 if( nelem == ipISO || dense.lgElmtOn[nelem] ) 00233 { 00234 for( i=0; i<iso.numLevels_max[ipISO][nelem]; i++ ) 00235 { 00236 Punch1LineData( &SatelliteLines[ipISO][nelem][i], ioPUN , true ); 00237 } 00238 } 00239 } 00240 } 00241 00242 /* want very small ionized fractions so get mostly H2 cs */ 00243 dense.eden = 1e-6; 00244 dense.gas_phase[ipHYDROGEN] = 1e-6f; 00245 dense.EdenHCorr = 1e-6f; 00246 dense.xIonDense[ipHYDROGEN][1] = 1.; 00247 hmi.Hmolec[ipMH2g] = 1.; 00248 hmi.Hmolec[ipMH2s] = 1.; 00249 dense.xIonDense[ipHYDROGEN][1] = 1e-6f; 00250 00251 /* H2 molecule */ 00252 fprintf( ioPUN, "\n\n\n" ); 00253 fprintf( ioPUN, " end database, start H2 lines\n" ); 00254 fprintf( ioPUN, "Eu Vu Ju El Vl Jl WL gl gu gf A CS n(crt)\n" ); 00255 00256 /* ioPUN unit, and option to print all possible lines - false indicates 00257 * punch only significant lines */ 00258 H2_LevelPops(); 00259 H2_Punch_line_data( ioPUN , false ); 00260 00261 fprintf( ioPUN, "\n\n\n" ); 00262 fprintf( ioPUN, " end H2, start 12CO rotation lines\n" ); 00263 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00264 CO_PopsEmisCool(&C12O16Rotate, nCORotate ,1. , "12CO",&a,&b ); 00265 for( j=0; j< nCORotate; ++j ) 00266 { 00267 Punch1LineData( &C12O16Rotate[j] , ioPUN , true); 00268 } 00269 00270 fprintf( ioPUN, "\n\n\n" ); 00271 fprintf( ioPUN, " end 12CO start 13CO rotation lines\n" ); 00272 fprintf( ioPUN, "Ion WL gl gu gf A CS n(crt)\n" ); 00273 CO_PopsEmisCool(&C13O16Rotate, nCORotate , 1. ,"13CO",&a,&b ); 00274 for( j=0; j< nCORotate; ++j ) 00275 { 00276 Punch1LineData( &C13O16Rotate[j] , ioPUN , true); 00277 } 00278 00279 /* punch FeII data if atom is currently enabled */ 00280 fprintf( ioPUN, "\n\n\n" ); 00281 fprintf( ioPUN, " end 13CO rotation lines, start FeII lines\n" ); 00282 fprintf( ioPUN, " Lo Hi Ion label WL gl gu gf A CS n(crt)\n" ); 00283 00284 /* ioPUN unit, and option to print all possible lines - false indicates 00285 * punch only significant lines */ 00286 FeIIPunData( ioPUN , false ); 00287 00288 /* stop when done, we have done some serious damage to the code */ 00289 cdEXIT(EXIT_SUCCESS); 00290 } 00291 00292 /*Punch1LineData punch data for one line */ 00293 void Punch1LineData( transition * t , FILE * ioPUN , 00294 /* flag saying whether to give collision strength too - in multi level atoms 00295 * it will be not valid without a great deal more work */ 00296 bool lgCS_2 ) 00297 { 00298 00299 double CritDen; 00300 /* these are used for parts of the line label */ 00301 char chLbl[11]; 00302 00303 DEBUG_ENTRY( "Punch1LineData()" ); 00304 00305 if( t->ipCont <= 0 ) 00306 { 00307 return; 00308 } 00309 00315 /*iWL = iWavLen( t , &chUnits , &chShift );*/ 00316 /* ion label, like C 1 */ 00317 chIonLbl(chLbl , t ); 00318 fprintf(ioPUN,"%s\t", chLbl ); 00319 00320 /* this is the second piece of the line label, pick up string after start */ 00321 00322 /* the wavelength */ 00323 prt_wl(ioPUN, t->WLAng ); 00324 00325 fprintf( ioPUN, " %3ld%3ld", 00326 /* lower and upper stat weights */ 00327 (long)(t->Lo->g), 00328 (long)(t->Hi->g) ); 00329 00330 /* oscillator strength */ 00331 fprintf( ioPUN,PrintEfmt("%9.2e", t->Emis->gf)); 00332 00333 /* Einstein A for transition */ 00334 fprintf( ioPUN,PrintEfmt("%9.2e", t->Emis->Aul)); 00335 00336 /* next collision strengths, use different formats depending on size 00337 * of collision strength */ 00338 if( t->Coll.cs > 100. ) 00339 { 00340 fprintf( ioPUN, "%7.1f", t->Coll.cs ); 00341 } 00342 else if( t->Coll.cs > 10. ) 00343 { 00344 fprintf( ioPUN, "%7.2f", t->Coll.cs ); 00345 } 00346 else if( t->Coll.cs > 1. ) 00347 { 00348 fprintf( ioPUN, "%7.3f", t->Coll.cs ); 00349 } 00350 else if( t->Coll.cs > .01 ) 00351 { 00352 fprintf( ioPUN, "%7.4f", t->Coll.cs ); 00353 } 00354 else if( t->Coll.cs > 0.0 ) 00355 { 00356 fprintf( ioPUN, " %.3e", t->Coll.cs ); 00357 } 00358 else 00359 { 00360 fprintf( ioPUN, "%7.4f", 0. ); 00361 } 00362 00363 /* now print critical density but only if cs is positive 00364 * >>chng 06 mar 24, add flag lgCS_2 - in multi-level systems do not want 00365 * to punch cs since not computed properly */ 00366 if( lgCS_2 && t->Coll.cs> 0. ) 00367 { 00368 CritDen = t->Emis->Aul * t->Hi->g*phycon.sqrte / (t->Coll.cs*COLL_CONST); 00369 CritDen = log10(CritDen); 00370 } 00371 else 00372 { 00373 CritDen = 0.; 00374 } 00375 fprintf( ioPUN, "%7.3f\n",CritDen ); 00376 return; 00377 }