cloudy: /home66/gary/public_html/cloudy/c13_branch/source/save

00001 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
00002  * others.  For conditions of distribution and use see copyright notice in license.txt */
00003 /*save_line parse save lines command, or actually do the save lines output */
00004 /*Save_Line_RT parse the save 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 "save.h"
00018 #include "iso.h"
00019 #include "parser.h"
00020 /* this is the limit to the number of emission lines we can store */
00021 #define NPUNLM  200L
00022 
00023 /* implement the save line xxx command.  cumulative, structure, and
00024  * emissivity all use same code base and variables, so only one can be used
00025  * at present */
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   /* true, return rel intensity, false, log of luminosity or intensity I */
00035   bool lgLog3,
00036   char *chHeader)
00037 {
00038         char chTemp[INPUT_LINE_LENGTH];
00039 
00040         // save return value of cdLine, 0 for success, -number of lines for fail
00041         long int i;
00042 
00043         DEBUG_ENTRY( "parse_save_line()" );
00044 
00045         /* very first time this routine is called, chDo is "READ" and we read
00046          * in lines from the input stream.  The line labels and wavelengths
00047          * are store locally, and output in later calls to this routine
00048          * following is flag saying whether to do relative intensity or
00049          * absolute emissivity */
00050         lgRelativeIntensity = lgLog3;
00051         
00052         /* number of lines we will save */
00053         nLinesEntered = 0;
00054         
00055         /* get the next line, and check for eof */
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         /* convert line to caps */
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                 /* this is total number stored so far */
00079                 ++nLinesEntered;
00080                 
00081                 /* get next line and check for eof */
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, /* the file we will write to */
00103   const char *chDo, 
00104   // intrinsic or emergent line emission?
00105   bool lgEmergent)
00106 {
00107         // save return value of cdLine, 0 for success, -number of lines for fail
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                 /* save lines emissivity command */
00121                 static bool lgMustGetLines=true,
00122                         lgBadLine=false;
00123                 lgBadLine = false;
00124                 static bool lgBadH2Line;
00125                 /* it is possible that we will get here after an initial temperature
00126                  * too high abort, and the line arrays will not have been defined.
00127                  * do no lines in this case.  must still do save so that there
00128                  * is not a missing line in the grid save output */
00129                 if( LineSave.nsum >0 )
00130                 {
00131                         lgBadH2Line = false;
00132                         lgBadLine = false;
00133                         /* save lines structure command */
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                                                 // missed line - ignore if H2 line since large model may be off
00142                                                 if( !h2.lgEnabled && strncmp( chPLab[i] , "H2  " , 4 )==0 )
00143                                                 {
00144                                                         static bool lgMustPrintFirstTime = true;
00145                                                         if( lgMustPrintFirstTime )
00146                                                         {
00147                                                                 /* it's an H2 line and H2 is not being done - ignore it */
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                                                         /* flag saying to ignore this line */
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                                 /* 0th line is dummy, can't be used, so this is safe */
00167                                 ASSERT( ipLine[i] > 0 || lgBadLine || lgBadH2Line || 
00168                                         /* this is case where we did not find line on previous iteration,
00169                                          * perhaps because H2 is not turned on, and -1 or -2 was
00170                                          * stored */
00171                                         (ipLine[i]<0&&!lgMustGetLines) );
00172                                 /* this version of cdEmis uses index, does not search, do not call if line could not be found */
00173                                 /* test on case where we abort before first zone is done
00174                                  * this happens in grid when temperature bounds of code
00175                                  * are exceeded.  In this case return small float */
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                 /* print depth */
00191                 fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
00192 
00193                 /* then print all line emissivity */
00194                 for( i=0; i < nLinesEntered; i++ )
00195                 {
00196                         /*lint -e644 dlum not initialized */
00197                         fprintf( ioPUN, "\t%.4f", log10( MAX2( SMALLFLOAT , dlum[i] ) ) );
00198                         /*lint +e644 dlum not initialized */
00199                 }
00200                 fprintf( ioPUN, "\n" );
00201         }
00202 
00203         else if( strcmp(chDo,"PUNC") == 0 )
00204         {
00205                 /* save lines cumulative command */
00206                 fprintf( ioPUN, "%.5e", radius.depth_mid_zone );
00207 
00208                 /* it is possible that we will get here after an initial temperature
00209                  * too high abort, and the line arrays will not have been defined.
00210                  * do no lines in this case.  must still do save so that there
00211                  * is not a missing line in the grid save output */
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                                         /* relative intensity case */
00220                                         a[i] = relint;
00221                                 else
00222                                         /* emissivity or luminosity case */
00223                                         a[i] = absint;
00224 
00225                                 if( nCdLineReturn<=0 )
00226                                 {
00227                                         /* missed line - ignore if H2 line */
00228                                         if( !h2.lgEnabled && strncmp( chPLab[i] , "H2  " , 4 )==0 )
00229                                         {
00230                                                 static bool lgMustPrintFirstTime = true;
00231                                                 if( lgMustPrintFirstTime )
00232                                                 {
00233                                                         /* it's an H2 line and H2 is not being done - ignore it */
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                                                 /* flag saying to ignore this line */
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 /*Save_Line_RT parse the save line rt command - read in a set of lines */
00291 void Parse_Save_Line_RT(Parser &p)
00292 {
00293         /* save line rt parameters */
00294         DEBUG_ENTRY( "Parse_Save_Line_RT()" );
00295 
00296         /* very first time this routine is called, chDo is "READ" and we read
00297          * in lines from the input stream.  The line labels and wavelengths
00298          * are store locally, and output in later calls to this routine */
00299         
00300         /* say that we must print the header */
00301         lgMustPrintHeader = true;
00302         
00303         /* get the next line, and check for eof */
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                 /* right now it just does lines in the iso sequences */
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                 /* now increment number of lines */
00337                 ++nLine;
00338                 
00339                 /* now get next line until we hit eof or the word END */
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         /* save line rt parameters */
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                         /* print info for header of file, line id and pump rate */
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                                 /* for now code only exists for H He like iso - assert this */
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                 /* index for line within continuum array */
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
/home66/gary/public_html/cloudy/c13_branch/source/save_line.cpp
