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