00001
00002
00003
00004
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 "save.h"
00018 #include "iso.h"
00019 #include "parser.h"
00020
00021 #define NPUNLM 200L
00022
00023
00024
00025
00026
00027 static char chPLab[NPUNLM][5];
00028 static long int nLinesEntered;
00029 static realnum wavelength[NPUNLM];
00030 static long int ipLine[NPUNLM];
00031 static bool lgRelativeIntensity;
00032
00033 void parse_save_line(Parser &p,
00034
00035 bool lgLog3,
00036 char *chHeader)
00037 {
00038 char chTemp[INPUT_LINE_LENGTH];
00039
00040
00041 long int i;
00042
00043 DEBUG_ENTRY( "parse_save_line()" );
00044
00045
00046
00047
00048
00049
00050 lgRelativeIntensity = lgLog3;
00051
00052
00053 nLinesEntered = 0;
00054
00055
00056 p.getline();
00057 if( p.m_lgEOF )
00058 {
00059 fprintf( ioQQQ,
00060 " Hit EOF while reading line list; use END to end list.\n" );
00061 cdEXIT(EXIT_FAILURE);
00062 }
00063
00064
00065
00066 while( p.strcmp("END") != 0 )
00067 {
00068 if( nLinesEntered >= NPUNLM )
00069 {
00070 fprintf( ioQQQ,
00071 " Too many lines have been entered; the limit is %ld. Increase variable NPUNLM in routine save_line.\n",
00072 NPUNLM );
00073 cdEXIT(EXIT_FAILURE);
00074 }
00075
00076
00077 strncpy( chPLab[nLinesEntered], p.getCommandRaw(4).c_str() , 4 );
00078
00079
00080 chPLab[nLinesEntered][4] = 0;
00081
00082
00083 wavelength[nLinesEntered] = (realnum)p.getWaveOpt();
00084
00085
00086 ++nLinesEntered;
00087
00088
00089 p.getline();
00090 if( p.m_lgEOF )
00091 {
00092 fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
00093 cdEXIT(EXIT_FAILURE);
00094 }
00095
00096 }
00097
00098 sprintf( chHeader, "depth");
00099 for( i=0; i < nLinesEntered; i++ )
00100 {
00101 sprintf( chTemp, "\t%s ", chPLab[i] );
00102 strcat( chHeader, chTemp );
00103 sprt_wl( chTemp, wavelength[i] );
00104 strcat( chHeader, chTemp );
00105 }
00106 strcat( chHeader, "\n" );
00107 }
00108
00109 void save_line(FILE * ioPUN,
00110 const char *chDo,
00111
00112 bool lgEmergent)
00113 {
00114
00115 long int nCdLineReturn;
00116 long int i;
00117 double a[32],
00118 absint,
00119 emiss,
00120 relint;
00121 double dlum[NPUNLM];
00122
00123 DEBUG_ENTRY( "save_line()" );
00124
00125 if( strcmp(chDo,"PUNS") == 0 )
00126 {
00127
00128 static bool lgMustGetLines=true,
00129 lgBadLine=false;
00130 lgBadLine = false;
00131 static bool lgBadH2Line;
00132
00133
00134
00135
00136 if( LineSave.nsum >0 )
00137 {
00138 lgBadH2Line = false;
00139 lgBadLine = false;
00140
00141 for( i=0; i < nLinesEntered; i++ )
00142 {
00143 if( nzone <= 1 && lgMustGetLines )
00144 {
00145 if( (ipLine[i] = cdEmis((char*)chPLab[i],wavelength[i],
00146 &emiss,lgEmergent)) <= 0 )
00147 {
00148
00149 if( !h2.lgH2ON && strncmp( chPLab[i] , "H2 " , 4 )==0 )
00150 {
00151 static bool lgMustPrintFirstTime = true;
00152 if( lgMustPrintFirstTime )
00153 {
00154
00155 fprintf( ioQQQ,"\nPROBLEM Did not find an H2 line, the large model is not "
00156 "included, so I will ignore it. Log intensity set to -30.\n" );
00157 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n\n");
00158 lgMustPrintFirstTime = false;
00159 }
00160
00161 ipLine[i] = -2;
00162 lgBadH2Line = true;
00163 }
00164 else
00165 {
00166 fprintf( ioQQQ, " save_line could not find line: %s %f\n",
00167 chPLab[i], wavelength[i] );
00168 ipLine[i] = -1;
00169 lgBadLine = true;
00170 }
00171 }
00172 }
00173
00174 ASSERT( ipLine[i] > 0 || lgBadLine || lgBadH2Line ||
00175
00176
00177
00178 (ipLine[i]<0&&!lgMustGetLines) );
00179
00180
00181
00182
00183 if( lgAbort && nzone <=1 )
00184 dlum[i] = SMALLFLOAT;
00185 else if( ipLine[i]>0 )
00186 cdEmis_ip(ipLine[i],&dlum[i],lgEmergent);
00187 else
00188 dlum[i] = 1e-30f;
00189 }
00190 if( lgBadLine )
00191 {
00192 cdEXIT(EXIT_FAILURE);
00193 }
00194 }
00195 lgMustGetLines = false;
00196
00197
00198 fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
00199
00200
00201 for( i=0; i < nLinesEntered; i++ )
00202 {
00203
00204 fprintf( ioPUN, "\t%.4f", log10( MAX2( SMALLFLOAT , dlum[i] ) ) );
00205
00206 }
00207 fprintf( ioPUN, "\n" );
00208 }
00209
00210 else if( strcmp(chDo,"PUNC") == 0 )
00211 {
00212
00213 fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
00214
00215
00216
00217
00218
00219 if( LineSave.nsum >0 )
00220 {
00221 for( i=0; i < nLinesEntered; i++ )
00222 {
00223 nCdLineReturn = cdLine((char*)chPLab[i],wavelength[i],
00224 &relint,&absint,lgEmergent);
00225 if( lgRelativeIntensity )
00226
00227 a[i] = relint;
00228 else
00229
00230 a[i] = absint;
00231
00232 if( nCdLineReturn<=0 )
00233 {
00234
00235 if( !h2.lgH2ON && strncmp( chPLab[i] , "H2 " , 4 )==0 )
00236 {
00237 static bool lgMustPrintFirstTime = true;
00238 if( lgMustPrintFirstTime )
00239 {
00240
00241 fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
00242 "included, so I will ignore it. Log intensity set to -30.\n" );
00243 fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
00244 lgMustPrintFirstTime = false;
00245 }
00246
00247 a[i] = -30.;
00248 absint = -30.;
00249 relint = -30.;
00250 }
00251 else
00252 {
00253 fprintf( ioQQQ, " save_line could not fine line: %s %f\n",
00254 chPLab[i], wavelength[i] );
00255 cdEXIT(EXIT_FAILURE);
00256 }
00257 }
00258 }
00259
00260 for( i=0; i < nLinesEntered; i++ )
00261 {
00262 fprintf( ioPUN, "\t%.4e", a[i] );
00263 }
00264 }
00265 fprintf( ioPUN, "\n" );
00266 }
00267
00268 else
00269 {
00270 fprintf( ioQQQ,
00271 " unrecognized key for save_line=%4.4s\n",
00272 chDo );
00273 cdEXIT(EXIT_FAILURE);
00274 }
00275 return;
00276 }
00277
00278 #define LIMLINE 10
00279 static long int line_RT_type[LIMLINE] =
00280 {LONG_MIN , LONG_MIN ,LONG_MIN , LONG_MIN ,LONG_MIN ,
00281 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
00282 line_RT_ipISO[LIMLINE] =
00283 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
00284 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
00285 line_RT_nelem[LIMLINE] =
00286 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
00287 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
00288 line_RT_ipHi[LIMLINE] =
00289 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
00290 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN },
00291 line_RT_ipLo[LIMLINE] =
00292 {LONG_MIN , LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,
00293 LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN ,LONG_MIN };
00294 static bool lgMustPrintHeader=true;
00295 static long int nLine=-1;
00296
00297
00298 void Parse_Save_Line_RT(Parser &p)
00299 {
00300
00301 DEBUG_ENTRY( "Parse_Save_Line_RT()" );
00302
00303
00304
00305
00306
00307
00308 lgMustPrintHeader = true;
00309
00310
00311 nLine = 0;
00312 p.getline();
00313 if( p.m_lgEOF )
00314 {
00315 fprintf( ioQQQ,
00316 " Hit EOF while reading line list; use END to end list.\n" );
00317 cdEXIT(EXIT_FAILURE);
00318 }
00319
00320 do
00321 {
00322 if(nLine>=LIMLINE )
00323 {
00324 fprintf(ioQQQ," PUNCH RT has too many lines - increase LIMLINE in save_line.cpp\n");
00325 cdEXIT(EXIT_FAILURE);
00326 }
00327
00328
00329 line_RT_type[nLine] = (long int)p.FFmtRead();
00330 line_RT_ipISO[nLine] = (long int)p.FFmtRead();
00331 line_RT_nelem[nLine] = (long int)p.FFmtRead();
00332 line_RT_ipHi[nLine] = (long int)p.FFmtRead();
00333 line_RT_ipLo[nLine] = (long int)p.FFmtRead();
00334
00335 if( p.lgEOL() )
00336 {
00337 fprintf( ioQQQ,
00338 " there must be five numbers on this line\n");
00339 p.PrintLine(ioQQQ);
00340 cdEXIT(EXIT_FAILURE);
00341 }
00342
00343
00344 ++nLine;
00345
00346
00347 p.getline();
00348 } while( !p.m_lgEOF && !p.nMatch( "END") );
00349 if( p.m_lgEOF )
00350 {
00351 fprintf( ioQQQ,
00352 " Save_Line_RT hit end of file looking for END of RT lines\n");
00353 p.PrintLine(ioQQQ);
00354 cdEXIT(EXIT_FAILURE);
00355 }
00356 }
00357
00358 void Save_Line_RT(
00359 FILE * ioPUN )
00360 {
00361
00362
00363 DEBUG_ENTRY( "Save_Line_RT()" );
00364
00365
00366 static char chLabel[LIMLINE][30];
00367 long int n;
00368 if( lgMustPrintHeader )
00369 {
00370 fprintf( ioPUN , "Line\tP(con,inc)\tAul\tgl\tgu\n");
00371 for( n=0; n<nLine; ++n )
00372 {
00373
00374 sprintf( chLabel[n] , "%s ",
00375 chLineLbl(&Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]]) );
00376 fprintf( ioPUN , "%s ", chLabel[n] );
00377 fprintf( ioPUN , "%.4e ",
00378 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->pump);
00379 fprintf( ioPUN , "%.4e ",
00380 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->Aul);
00381 fprintf( ioPUN , "%.0f ",
00382 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Lo->g);
00383 fprintf( ioPUN , "%.0f ",
00384 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Hi->g);
00385 fprintf( ioPUN , "\n");
00386
00387 if( line_RT_type[n]!=0. )
00388 {
00389
00390 fprintf( ioQQQ,
00391 " PunchLine_RT only H, He like allowed for now\n");
00392 cdEXIT(EXIT_FAILURE);
00393 }
00394 }
00395 fprintf( ioPUN , "Line\tTauIn\tPopLo\tPopHi\tCul\tk(line)\tk(con,abs)\tk(con,scat)\n");
00396 lgMustPrintHeader = false;
00397 }
00398
00399 fprintf(ioPUN, "RADIUS\t%e\tDEPTH\t%e\tTe\t%e\tNe\t%e\n",
00400 radius.Radius_mid_zone ,
00401 radius.depth_mid_zone ,
00402 phycon.te ,
00403 dense.eden );
00404 for( n=0; n<nLine; ++n )
00405 {
00406
00407 long int ipCont = Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].ipCont;
00408 fprintf( ioPUN , "%s ", chLabel[n] );
00409 fprintf( ioPUN , "\t%e\t%e\t%e",
00410 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->TauIn ,
00411 StatesElemNEW[line_RT_nelem[n]][line_RT_nelem[n]-line_RT_ipISO[n]][line_RT_ipLo[n]].Pop ,
00412 StatesElemNEW[line_RT_nelem[n]][line_RT_nelem[n]-line_RT_ipISO[n]][line_RT_ipHi[n]].Pop
00413 );
00414 fprintf( ioPUN , "\t%e",
00415 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Coll.ColUL
00416 );
00417
00418 fprintf( ioPUN , "\t%e\t%e\t%e\n",
00419 Transitions[line_RT_ipISO[n]][line_RT_nelem[n]][line_RT_ipHi[n]][line_RT_ipLo[n]].Emis->PopOpc ,
00420 opac.opacity_abs[ipCont-1] ,
00421 opac.opacity_sct[ipCont-1]
00422 );
00423 }
00424 }
00425
00426 # undef LIMELM
00427