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