cloudy: /home66/gary/public_html/cloudy/c08_branch/source/punch

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
/home66/gary/public_html/cloudy/c08_branch/source/punch_line.cpp
