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 /*punch_line parse punch lines command, or actually do the punch lines output */ 00004 /*Punch_Line_RT parse the punch line rt command - read in a set of lines */ 00005 #include "cddefines.h" 00006 #include "cddrive.h" 00007 #include "radius.h" 00008 #include "taulines.h" 00009 #include "opacity.h" 00010 #include "phycon.h" 00011 #include "dense.h" 00012 #include "lines.h" 00013 #include "h2.h" 00014 #include "lines_service.h" 00015 #include "input.h" 00016 #include "prt.h" 00017 #include "punch.h" 00018 #include "iso.h" 00019 /* this is the limit to the number of emission lines we can store */ 00020 #define NPUNLM 200L 00021 00022 /* implement the punch line xxx command. cumulative, structure, and 00023 * emissivity all use same code base and variables, so only one can be used 00024 * at present */ 00025 void punch_line(FILE * ioPUN, /* the file we will write to */ 00026 const char *chDo, 00027 /* true, return rel intensity, false, log of luminosity or intensity I */ 00028 bool lgLog3, 00029 char *chHeader) 00030 { 00031 char chCap[INPUT_LINE_LENGTH], 00032 chCard[INPUT_LINE_LENGTH], 00033 chTemp[INPUT_LINE_LENGTH]; 00034 00035 static char chPLab[NPUNLM][5]; 00036 static long int nLinesEntered; 00037 static realnum wavelength[NPUNLM]; 00038 static long int ipLine[NPUNLM]; 00039 00040 bool lgEOF, 00041 lgEOL; 00042 // save return value of cdLine, 0 for success, -number of lines for fail 00043 long int nCdLineReturn; 00044 static bool lgRelativeIntensity; 00045 long int i; 00046 double a[32], 00047 absint, 00048 emiss, 00049 relint; 00050 double dlum[NPUNLM]; 00051 00052 DEBUG_ENTRY( "punch_line()" ); 00053 00054 if( strcmp(chDo,"READ") == 0 ) 00055 { 00056 /* very first time this routine is called, chDo is "READ" and we read 00057 * in lines from the input stream. The line labels and wavelengths 00058 * are store locally, and output in later calls to this routine 00059 * following is flag saying whether to do relative intensity or 00060 * absolute emissivity */ 00061 lgRelativeIntensity = lgLog3; 00062 00063 /* number of lines we will save */ 00064 nLinesEntered = 0; 00065 00066 /* get the next line, and check for eof */ 00067 input_readarray(chCard,&lgEOF); 00068 if( lgEOF ) 00069 { 00070 fprintf( ioQQQ, 00071 " Hit EOF while reading line list; use END to end list.\n" ); 00072 cdEXIT(EXIT_FAILURE); 00073 } 00074 00075 /* convert line to caps */ 00076 strcpy( chCap, chCard ); 00077 caps(chCap); 00078 00079 while( strncmp(chCap, "END" ,3 ) != 0 ) 00080 { 00081 if( nLinesEntered >= NPUNLM ) 00082 { 00083 fprintf( ioQQQ, 00084 " Too many lines have been entered; the limit is %ld. Increase variable NPUNLM in routine punch_line.\n", 00085 NPUNLM ); 00086 cdEXIT(EXIT_FAILURE); 00087 } 00088 00089 /* order on line is label (col 1-4), wavelength */ 00090 strncpy( chPLab[nLinesEntered], chCard , 4 ); 00091 00092 /* null terminate the string*/ 00093 chPLab[nLinesEntered][4] = 0; 00094 00095 /* now get wavelength */ 00096 i = 5; 00097 wavelength[nLinesEntered] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00098 00099 /* now convert wavelength to angstroms */ 00100 /* line was entered, look for possible micron or cm label */ 00101 if( input.chCARDCAPS[i-1] == 'M' ) 00102 { 00103 /* microns */ 00104 wavelength[nLinesEntered] *= 1e4; 00105 } 00106 else if( input.chCARDCAPS[i-1] == 'C' ) 00107 { 00108 /* microns */ 00109 wavelength[nLinesEntered] *= 1e8; 00110 } 00111 00112 /* this is total number stored so far */ 00113 ++nLinesEntered; 00114 00115 /* get next line and check for eof */ 00116 input_readarray(chCard,&lgEOF); 00117 if( lgEOF ) 00118 { 00119 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" ); 00120 cdEXIT(EXIT_FAILURE); 00121 } 00122 00123 /* convert it to caps */ 00124 strcpy( chCap, chCard ); 00125 caps(chCap); 00126 } 00127 00128 /*fprintf( ioPUN, "%li lines were entered, they were;\n", 00129 nLinesEntered );*/ 00130 00131 /* >>chng 06 nov 17, this check needed so headers are only printed once in grid mode. */ 00132 sprintf( chHeader, "#depth\t"); 00133 for( i=0; i < nLinesEntered; i++ ) 00134 { 00135 sprintf( chTemp, "%s ", chPLab[i] ); 00136 strcat( chHeader, chTemp ); 00137 sprt_wl( chTemp, wavelength[i] ); 00138 strcat( chHeader, chTemp ); 00139 strcat( chHeader, "\t" ); 00140 } 00141 strcat( chHeader, "\n" ); 00142 } 00143 00144 else if( strcmp(chDo,"PUNS") == 0 ) 00145 { 00146 /* punch lines emissivity command */ 00147 static bool lgMustGetLines=true, 00148 lgBadLine=false; 00149 lgBadLine = false; 00150 static bool lgBadH2Line; 00151 /* it is possible that we will get here after an initial temperature 00152 * too high abort, and the line arrays will not have been defined. 00153 * do no lines in this case. must still do punch so that there 00154 * is not a missing line in the grid punch output */ 00155 if( LineSave.nsum >0 ) 00156 { 00157 lgBadH2Line = false; 00158 lgBadLine = false; 00159 /* punch lines structure command */ 00160 for( i=0; i < nLinesEntered; i++ ) 00161 { 00162 if( nzone <= 1 && lgMustGetLines ) 00163 { 00164 if( (ipLine[i] = cdEmis((char*)chPLab[i],wavelength[i],&emiss)) <= 0 ) 00165 { 00166 /* missed line - ignore if H2 line */ 00167 if( !h2.lgH2ON && strncmp( chPLab[i] , "H2 " , 4 )==0 ) 00168 { 00169 static bool lgMustPrintFirstTime = true; 00170 if( lgMustPrintFirstTime ) 00171 { 00172 /* it's an H2 line and H2 is not being done - ignore it */ 00173 fprintf( ioQQQ,"\nPROBLEM Did not find an H2 line, the large model is not " 00174 "included, so I will ignore it. Log intensity set to -30.\n" ); 00175 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n\n"); 00176 lgMustPrintFirstTime = false; 00177 } 00178 /* flag saying to ignore this line */ 00179 ipLine[i] = -1; 00180 lgBadH2Line = true; 00181 } 00182 else 00183 { 00184 fprintf( ioQQQ, " PUNLIN could not find line: %s %f\n", 00185 chPLab[i], wavelength[i] ); 00186 ipLine[i] = -1; 00187 lgBadLine = true; 00188 } 00189 } 00190 } 00191 /* 0th line is dummy, can't be used, so this is safe */ 00192 ASSERT( ipLine[i] > 0 || lgBadLine || lgBadH2Line ); 00193 /* this version of cdEmis uses index, does not search, do not call if line could not be found */ 00194 /* test on case where we abort before first zone is done 00195 * this happens in grid when temperature bounds of code 00196 * are exceeded. In this case return small float */ 00197 if( lgAbort && nzone <=1 ) 00198 dlum[i] = SMALLFLOAT; 00199 else if( ipLine[i]>0 ) 00200 cdEmis_ip(ipLine[i],&dlum[i]); 00201 else 00202 dlum[i] = 1e-30f; 00203 } 00204 if( lgBadLine ) 00205 { 00206 cdEXIT(EXIT_FAILURE); 00207 } 00208 } 00209 lgMustGetLines = false; 00210 00211 /* print depth */ 00212 fprintf( ioPUN, "%.5e", radius.depth_mid_zone ); 00213 00214 /* then print all line emissivity */ 00215 for( i=0; i < nLinesEntered; i++ ) 00216 { 00217 /*lint -e644 dlum not initialized */ 00218 fprintf( ioPUN, "\t%.4f", log10( MAX2( SMALLFLOAT , dlum[i] ) ) ); 00219 /*lint +e644 dlum not initialized */ 00220 } 00221 fprintf( ioPUN, "\n" ); 00222 } 00223 00224 else if( strcmp(chDo,"PUNC") == 0 ) 00225 { 00226 /* punch lines cumulative command */ 00227 fprintf( ioPUN, "%.5e", radius.depth_mid_zone ); 00228 00229 /* it is possible that we will get here after an initial temperature 00230 * too high abort, and the line arrays will not have been defined. 00231 * do no lines in this case. must still do punch so that there 00232 * is not a missing line in the grid punch output */ 00233 if( LineSave.nsum >0 ) 00234 { 00235 for( i=0; i < nLinesEntered; i++ ) 00236 { 00237 if( lgRelativeIntensity ) 00238 { 00239 /* relative intensity case */ 00240 nCdLineReturn = cdLine((char*)chPLab[i],wavelength[i],&a[i],&absint); 00241 } 00242 else 00243 { 00244 /* emissivity or luminosity case */ 00245 nCdLineReturn = cdLine((char*)chPLab[i],wavelength[i],&relint,&a[i]); 00246 } 00247 00248 if( nCdLineReturn<=0 ) 00249 { 00250 /* missed line - ignore if H2 line */ 00251 if( !h2.lgH2ON && strncmp( chPLab[i] , "H2 " , 4 )==0 ) 00252 { 00253 static bool lgMustPrintFirstTime = true; 00254 if( lgMustPrintFirstTime ) 00255 { 00256 /* it's an H2 line and H2 is not being done - ignore it */ 00257 fprintf( ioQQQ,"Did not find an H2 line, the large model is not " 00258 "included, so I will ignore it. Log intensity set to -30.\n" ); 00259 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n"); 00260 lgMustPrintFirstTime = false; 00261 } 00262 /* flag saying to ignore this line */ 00263 a[i] = -30.; 00264 absint = -30.; 00265 relint = -30.; 00266 } 00267 else 00268 { 00269 fprintf( ioQQQ, " PUNLIN could not fine line: %s %f\n", 00270 chPLab[i], wavelength[i] ); 00271 cdEXIT(EXIT_FAILURE); 00272 } 00273 } 00274 } 00275 00276 for( i=0; i < nLinesEntered; i++ ) 00277 { 00278 fprintf( ioPUN, "\t%.4e", a[i] ); 00279 } 00280 } 00281 fprintf( ioPUN, "\n" ); 00282 } 00283 00284 else 00285 { 00286 fprintf( ioQQQ, 00287 " unrecognized key for punch_line=%4.4s\n", 00288 chDo ); 00289 cdEXIT(EXIT_FAILURE); 00290 } 00291 return; 00292 } 00293 00294 /*Punch_Line_RT parse the punch line rt command - read in a set of lines */ 00295 void Punch_Line_RT( 00296 FILE * ioPUN, /* the file we will write to */ 00297 const char *chDo ) 00298 { 00299 # define LIMLINE 10 00300 /* punch line rt parameters */ 00301 static long int 00302 line_RT_type[LIMLINE] = 00303 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN }, 00304 line_RT_ipISO[LIMLINE] = 00305 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN }, 00306 line_RT_nelem[LIMLINE] = 00307 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN }, 00308 line_RT_ipHi[LIMLINE] = 00309 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN }, 00310 line_RT_ipLo[LIMLINE] = 00311 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN }; 00312 static bool lgMustPrintHeader=true; 00313 static long int nLine=-1; 00314 long int i; 00315 bool lgEOF , lgEOL; 00316 char chCap[LIMLINE][INPUT_LINE_LENGTH], 00317 chCard[INPUT_LINE_LENGTH]; 00318 00319 00320 DEBUG_ENTRY( "Punch_Line_RT()" ); 00321 00322 if( strcmp(chDo,"READ") == 0 ) 00323 { 00324 /* very first time this routine is called, chDo is "READ" and we read 00325 * in lines from the input stream. The line labels and wavelengths 00326 * are store locally, and output in later calls to this routine */ 00327 00328 /* say that we must print the header */ 00329 lgMustPrintHeader = true; 00330 00331 /* get the next line, and check for eof */ 00332 nLine = 0; 00333 input_readarray(chCard,&lgEOF); 00334 if( lgEOF ) 00335 { 00336 fprintf( ioQQQ, 00337 " Hit EOF while reading line list; use END to end list.\n" ); 00338 cdEXIT(EXIT_FAILURE); 00339 } 00340 00341 do 00342 { 00343 if(nLine>=LIMLINE ) 00344 { 00345 fprintf(ioQQQ," PUNCH RT has too many lines - increase LIMLINE in punch_line.c\n"); 00346 cdEXIT(EXIT_FAILURE); 00347 } 00348 00349 /* convert line to caps */ 00350 strcpy( chCap[nLine], chCard ); 00351 caps(chCap[nLine]); 00352 00353 /* right now it just does lines in the iso sequences */ 00354 i = 1; 00355 line_RT_type[nLine] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00356 line_RT_ipISO[nLine] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00357 line_RT_nelem[nLine] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00358 line_RT_ipHi[nLine] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00359 line_RT_ipLo[nLine] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00360 00361 if( lgEOL ) 00362 { 00363 fprintf( ioQQQ, 00364 " there must be five numbers on this line =%s\n", 00365 chCard ); 00366 cdEXIT(EXIT_FAILURE); 00367 } 00368 00369 /* now increment number of lines */ 00370 ++nLine; 00371 00372 /* now get next line until we hit eof or the word END */ 00373 input_readarray(chCard,&lgEOF); 00374 caps(chCard); 00375 } while( !lgEOF && !nMatch( "END" , chCard) ); 00376 if( lgEOF ) 00377 { 00378 fprintf( ioQQQ, 00379 " Punch_Line_RT hit end of file looking for END of RT lines =%s\n", 00380 chCard ); 00381 cdEXIT(EXIT_FAILURE); 00382 } 00383 } 00384 00385 else if( strcmp(chDo,"PUNC") == 0 ) 00386 { 00387 static char chLabel[LIMLINE][30]; 00388 long int n; 00389 if( lgMustPrintHeader ) 00390 { 00391 fprintf( ioPUN , "Line\tP(con,inc)\tAul\tgl\tgu\n"); 00392 for( n=0; n<nLine; ++n ) 00393 { 00394 /* print info for header of file, line id and pump rate */ 00395 sprintf( chLabel[n] , "%s ", 00396 chLineLbl(&Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]]) ); 00397 fprintf( ioPUN , "%s ", chLabel[n] ); 00398 fprintf( ioPUN , "%.4e ", 00399 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->pump); 00400 fprintf( ioPUN , "%.4e ", 00401 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->Aul); 00402 fprintf( ioPUN , "%.0f ", 00403 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Lo->g); 00404 fprintf( ioPUN , "%.0f ", 00405 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Hi->g); 00406 fprintf( ioPUN , "\n"); 00407 00408 if( line_RT_type[n]!=0. ) 00409 { 00410 /* for now code only exists for H He like iso - assert this */ 00411 fprintf( ioQQQ, 00412 " PunchLine_RT only H, He like allowed for now\n"); 00413 cdEXIT(EXIT_FAILURE); 00414 } 00415 } 00416 fprintf( ioPUN , "Line\tTauIn\tPopLo\tPopHi\tCul\tk(line)\tk(con,abs)\tk(con,scat)\n"); 00417 lgMustPrintHeader = false; 00418 } 00419 00420 fprintf(ioPUN, "RADIUS\t%e\tDEPTH\t%e\tTe\t%e\tNe\t%e\n", 00421 radius.Radius_mid_zone , 00422 radius.depth_mid_zone , 00423 phycon.te , 00424 dense.eden ); 00425 for( n=0; n<nLine; ++n ) 00426 { 00427 /* index for line within continuum array */ 00428 long int ipCont = Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].ipCont; 00429 fprintf( ioPUN , "%s ", chLabel[n] ); 00430 fprintf( ioPUN , "\t%e\t%e\t%e", 00431 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->TauIn , 00432 StatesElem[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipLo[n]].Pop , 00433 StatesElem[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]].Pop 00434 ); 00435 fprintf( ioPUN , "\t%e", 00436 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Coll.ColUL 00437 ); 00438 00439 fprintf( ioPUN , "\t%e\t%e\t%e\n", 00440 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->PopOpc , 00441 opac.opacity_abs[ipCont-1] , 00442 opac.opacity_sct[ipCont-1] 00443 ); 00444 } 00445 } 00446 00447 else 00448 { 00449 fprintf( ioQQQ, 00450 " internal error, unrecognized key for punch line rt=%4.4s\n", 00451 chDo ); 00452 cdEXIT(EXIT_FAILURE); 00453 } 00454 00455 return; 00456 } 00457 # undef LIMELM 00458