cloudy trunk
|
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 /*cdDrive main routine to call cloudy under all circumstances) */ 00004 /*cdReasonGeo write why the model stopped and type of geometry on io file */ 00005 /*cdWarnings write all warnings entered into comment stack */ 00006 /*cdEms obtain the local emissivity for a line, for the last computed zone */ 00007 /*cdColm get the column density for a constituent */ 00008 /*cdLine get the predicted line intensity, also index for line in stack */ 00009 /*cdLine_ip get the predicted line intensity, using index for line in stack */ 00010 /*cdDLine get the predicted emergent line intensity, also index for line in stack */ 00011 /*cdCautions print out all cautions after calculation, on arbitrary io unit */ 00012 /*cdTemp_last routine to query results and return temperature of last zone */ 00013 /*cdDepth_depth get depth structure from previous iteration */ 00014 /*cdTimescales returns thermal, recombination, and H2 formation timescales */ 00015 /*cdSurprises print out all surprises on arbitrary unit number */ 00016 /*cdNotes print stack of notes about current calculation */ 00017 /*cdPressure_last routine to query results and return pressure of last zone */ 00018 /*cdTalk tells the code whether to print results or be silent */ 00019 /*cdOutp redirect output to arbitrary Fortran unit number */ 00020 /*cdRead routine to read in command lines when cloudy used as subroutine */ 00021 /*cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit */ 00022 /*cdIonFrac get ionization fractions for a constituent */ 00023 /*cdTemp get mean electron temperature for any element */ 00024 /*cdCooling_last routine to query results and return cooling of last zone */ 00025 /*cdHeating_last routine to query results and return heating of last zone */ 00026 /*cdEDEN_last return electron density of last zone */ 00027 /*cdNoExec call this routine to tell code not to actually execute */ 00028 /*cdDate - puts date of code into string */ 00029 /*cdVersion produces string that gives version number of the code */ 00030 /*cdExecTime any routine can call this, find the time [s] since cdInit was called */ 00031 /*cdPrintCommands( FILE *) prints all input commands into file */ 00032 /*cdDrive main routine to call cloudy under all circumstances) */ 00033 /*cdNwcns get the number of cautions and warnings, to tell if calculation is ok */ 00034 /*debugLine provides a debugging hook into the main line array */ 00035 /*cdEms_ip obtain the local emissivity for a line with known index */ 00036 /*cdnZone gets number of zones */ 00037 /*cdClosePunchFiles closes all the punch files that have been used */ 00038 /*cdLineListPunch create a file with a list of all emission lines punched, 00039 *and their index within the emission line stack */ 00040 /*cdB21cm - returns B as measured by 21 cm */ 00041 /*cdPrtWL print line wavelengths in Angstroms in the standard format */ 00042 #include "cddefines.h" 00043 #include "trace.h" 00044 #include "conv.h" 00045 #include "init.h" 00046 #include "lines.h" 00047 #include "pressure.h" 00048 #include "prt.h" 00049 #include "colden.h" 00050 #include "dense.h" 00051 #include "radius.h" 00052 #include "struc.h" 00053 #include "mole.h" 00054 #include "elementnames.h" 00055 #include "mean.h" 00056 #include "phycon.h" 00057 #include "called.h" 00058 #include "parse.h" 00059 #include "input.h" 00060 #include "taulines.h" 00061 #include "version.h" 00062 #include "thermal.h" 00063 #include "optimize.h" 00064 #include "grid.h" 00065 #include "timesc.h" 00066 #include "cloudy.h" 00067 #include "warnings.h" 00068 #include "lines_service.h" 00069 #include "cddrive.h" 00070 00071 /************************************************************************* 00072 * 00073 * cdDrive - main routine to call cloudy - returns 0 if all ok, 1 if problems 00074 * 00075 ************************************************************************/ 00076 00077 int cdDrive(void ) 00078 { 00079 bool lgBAD; 00080 00081 DEBUG_ENTRY( "cdDrive()" ); 00082 /********************************* 00083 * main routine to call cloudy * 00084 *********************************/ 00085 00086 /* this is set false when code loaded, set true when cdInit called, 00087 * this is check that cdInit was called first */ 00088 if( !lgcdInitCalled ) 00089 { 00090 printf(" cdInit was not called first - this must be the first call.\n"); 00091 cdEXIT(EXIT_FAILURE); 00092 } 00093 00094 if( trace.lgTrace ) 00095 { 00096 fprintf( ioQQQ, 00097 "cdDrive: lgOptimr=%1i lgVaryOn=%1i lgNoVary=%1i input.nSave:%li\n", 00098 optimize.lgOptimr , optimize.lgVaryOn , optimize.lgNoVary, input.nSave ); 00099 } 00100 00101 /* should we call cloudy, or the optimization driver? 00102 * possible to have VARY on line without OPTIMIZE being set 00103 * lgNoVary set with "no optimize" command */ 00104 if( optimize.lgOptimr && optimize.lgVaryOn && !optimize.lgNoVary ) 00105 optimize.lgVaryOn = true; 00106 else 00107 optimize.lgVaryOn = false; 00108 00109 /* one time initialization of core load - returns if already called 00110 * called here rather than in cdInit since at this point we know if 00111 * single sim or grid */ 00112 InitCoreload(); 00113 /* count now many sims have been done in this coreload, above set to 00114 * zero if first call, does nothing on subsequent calls 00115 * this increment brings to 1 on first sim */ 00116 ++nSimThisCoreload; 00117 00118 if( optimize.lgVaryOn ) 00119 { 00120 /* this branch called if optimizing or grid calculation */ 00121 if( trace.lgTrace ) 00122 fprintf( ioQQQ, "cdDrive: calling grid_do\n"); 00123 /* option to drive optimizer set if OPTIMIZE was in input stream */ 00124 lgBAD = grid_do(); 00125 } 00126 else 00127 { 00128 if( trace.lgTrace ) 00129 fprintf( ioQQQ, "cdDrive: calling cloudy\n"); 00130 try { 00131 00132 /* optimize did not occur, only compute one model, call cloudy */ 00133 lgBAD = cloudy(); 00134 } 00135 catch( bad_mpi &bad_info ) 00136 { 00137 fprintf( ioQQQ, "An MPI run failed. The fail mode is %li\n", bad_info.failMode() ); 00138 ClosePunchFiles( false ); 00139 lgAbort = true; 00140 lgBAD = true; 00141 } 00142 } 00143 00144 /* reset flag saying that cdInit has not been called */ 00145 lgcdInitCalled = false; 00146 00147 if( lgAbort || lgBAD ) 00148 { 00149 if( trace.lgTrace ) 00150 fprintf( ioQQQ, "cdDrive: returning failure during call. \n"); 00151 /* lgAbort set true if something wrong, so return lgBAD false. */ 00152 return(1); 00153 } 00154 else 00155 { 00156 /* everything is ok, return 0 */ 00157 return(0); 00158 } 00159 } 00160 00161 /************************************************************************* 00162 * 00163 * cdPrtWL write emission line wavelength 00164 * 00165 ************************************************************************/ 00166 00167 /*cdPrtWL print line wavelengths in Angstroms in the standard format - 00168 * a wrapper for prt_wl which must be kept parallel with sprt_wl 00169 * both of those live in pdt.c */ 00170 void cdPrtWL( FILE *io , realnum wl ) 00171 { 00172 00173 DEBUG_ENTRY( "cdPrtWL()" ); 00174 00175 prt_wl( io , wl ); 00176 return; 00177 } 00178 00179 00180 /************************************************************************* 00181 * 00182 * cdReasonGeo write why the model stopped and type of geometry on io file 00183 * 00184 ************************************************************************/ 00185 00186 00187 /*cdReasonGeo write why the model stopped and type of geometry on io file */ 00188 void cdReasonGeo(FILE * ioOUT) 00189 { 00190 00191 DEBUG_ENTRY( "cdReasonGeo()" ); 00192 00193 /*this is the reason the calculation stopped*/ 00194 fprintf( ioOUT, "%s", warnings.chRgcln[0] ); 00195 fprintf( ioOUT , "\n" ); 00196 /* this is the geometry */ 00197 fprintf( ioOUT, "%s", warnings.chRgcln[1] ); 00198 fprintf( ioOUT , "\n" ); 00199 return; 00200 } 00201 00202 00203 /************************************************************************* 00204 * 00205 * cdWarnings write all warnings entered into comment stack 00206 * 00207 ************************************************************************/ 00208 00209 /*cdWarnings write all warnings entered into comment stack */ 00210 00211 void cdWarnings(FILE *ioPNT ) 00212 { 00213 long int i; 00214 00215 DEBUG_ENTRY( "cdWarnings()" ); 00216 00217 if( warnings.nwarn > 0 ) 00218 { 00219 for( i=0; i < warnings.nwarn; i++ ) 00220 { 00221 /* these are all warnings that were entered in comment */ 00222 fprintf( ioPNT, "%s", warnings.chWarnln[i] ); 00223 fprintf( ioPNT, "\n" ); 00224 } 00225 } 00226 00227 return; 00228 } 00229 00230 00231 /************************************************************************* 00232 * 00233 * cdCautions print out all cautions after calculation, on arbitrary io unit 00234 * 00235 ************************************************************************/ 00236 00237 /*cdCautions print out all cautions after calculation, on arbitrary io unit */ 00238 00239 void cdCautions(FILE * ioOUT) 00240 { 00241 long int i; 00242 00243 DEBUG_ENTRY( "cdCautions()" ); 00244 00245 if( warnings.ncaun > 0 ) 00246 { 00247 for( i=0; i < warnings.ncaun; i++ ) 00248 { 00249 fprintf( ioOUT, "%s", warnings.chCaunln[i] ); 00250 fprintf( ioOUT, "\n" ); 00251 } 00252 } 00253 return; 00254 } 00255 00256 /************************************************************************* 00257 * 00258 * cdTimescales returns thermal, recombination, and H2 formation timescales 00259 * 00260 ************************************************************************/ 00261 00262 void cdTimescales( 00263 /* the thermal cooling timescale */ 00264 double *TTherm , 00265 /* the hydrogen recombination timescale */ 00266 double *THRecom , 00267 /* the H2 formation timescale */ 00268 double *TH2 ) 00269 { 00270 00271 DEBUG_ENTRY( "cdTimescales()" ); 00272 00273 /* these were all evaluated in AgeCheck, which was called by PrtComment */ 00274 00275 /* thermal or cooling timescale */ 00276 *TTherm = timesc.time_therm_long; 00277 00278 /* the hydrogen recombination timescale */ 00279 *THRecom = timesc.time_Hrecom_long; 00280 00281 /* longer of the the H2 formation and destruction timescales */ 00282 *TH2 = MAX2( timesc.time_H2_Dest_longest , timesc.time_H2_Form_longest ); 00283 return; 00284 } 00285 00286 00287 /************************************************************************* 00288 * 00289 * cdSurprises print out all surprises on arbitrary unit number 00290 * 00291 ************************************************************************/ 00292 00293 /*cdSurprises print out all surprises on arbitrary unit number */ 00294 00295 void cdSurprises(FILE * ioOUT) 00296 { 00297 long int i; 00298 00299 DEBUG_ENTRY( "cdSurprises()" ); 00300 00301 if( warnings.nbang > 0 ) 00302 { 00303 for( i=0; i < warnings.nbang; i++ ) 00304 { 00305 fprintf( ioOUT, "%s", warnings.chBangln[i] ); 00306 fprintf( ioOUT, "\n" ); 00307 } 00308 } 00309 00310 return; 00311 } 00312 00313 00314 /************************************************************************* 00315 * 00316 * cdNotes print stack of notes about current calculation 00317 * 00318 ************************************************************************/ 00319 00320 /*cdNotes print stack of notes about current calculation */ 00321 00322 void cdNotes(FILE * ioOUT) 00323 { 00324 long int i; 00325 00326 DEBUG_ENTRY( "cdNotes()" ); 00327 00328 if( warnings.nnote > 0 ) 00329 { 00330 for( i=0; i < warnings.nnote; i++ ) 00331 { 00332 fprintf( ioOUT, "%s", warnings.chNoteln[i] ); 00333 fprintf( ioOUT, "\n" ); 00334 } 00335 } 00336 return; 00337 } 00338 00339 /************************************************************************* 00340 * 00341 * cdCooling_last routine to query results and return cooling of last zone 00342 * 00343 ************************************************************************/ 00344 00345 /*cdCooling_last routine to query results and return cooling of last zone */ 00346 double cdCooling_last(void) /* return cooling for last zone */ 00347 { 00348 return (thermal.ctot); 00349 } 00350 00351 /************************************************************************* 00352 * 00353 * cdVersion - puts version number of code into string 00354 * incoming string must have sufficient length and will become null 00355 * terminated string 00356 * 00357 ************************************************************************/ 00358 00359 void cdVersion(char chString[]) 00360 { 00361 strcpy( chString , t_version::Inst().chVersion ); 00362 return; 00363 } 00364 00365 /************************************************************************* 00366 * 00367 * cdDate - puts date of code into string 00368 * incoming string must have at least 8 char and will become null 00369 * terminated string 00370 * 00371 ************************************************************************/ 00372 00373 /* cdDate - puts date of code into string */ 00374 void cdDate(char chString[]) 00375 { 00376 strcpy( chString , t_version::Inst().chDate ); 00377 return; 00378 } 00379 00380 00381 /************************************************************************* 00382 * 00383 * cdHeating_last routine to query results and return heating of last zone 00384 * 00385 ************************************************************************/ 00386 00387 /*cdHeating_last routine to query results and return heating of last zone */ 00388 00389 double cdHeating_last(void) /* return heating for last zone */ 00390 { 00391 return (thermal.htot); 00392 } 00393 00394 00395 /************************************************************************* 00396 * 00397 * cdEDEN_last return electron density of last zone 00398 * 00399 ************************************************************************/ 00400 00401 /*cdEDEN_last return electron density of last zone */ 00402 00403 double cdEDEN_last(void) /* return electron density for last zone */ 00404 { 00405 return ( dense.eden ); 00406 } 00407 00408 /************************************************************************* 00409 * 00410 * cdNoExec call this routine to tell code not to actually execute 00411 * 00412 ************************************************************************/ 00413 00414 /*cdNoExec call this routine to tell code not to actually execute */ 00415 #include "noexec.h" 00416 00417 void cdNoExec(void) 00418 { 00419 00420 DEBUG_ENTRY( "cdNoExec()" ); 00421 00422 /* option to read in all input quantiles and NOT execute the actual model 00423 * only check on input parameters - set by calling cdNoExec */ 00424 noexec.lgNoExec = true; 00425 00426 return; 00427 } 00428 00429 00430 /************************************************************************* 00431 * 00432 * cdSetExecTime routine to initialize variables keeping track of time at start of calculation 00433 * 00434 ************************************************************************/ 00435 00436 /* set to false initially, then to true when cdSetExecTime() is called to 00437 * initialize the clock */ 00438 static bool lgCalled=false; 00439 00440 /* >>chng 06 dec 19, RP rm "|| defined(__HP_aCC)" to run on hp */ 00441 #if defined(_MSC_VER) 00442 /* _MSC_VER branches assume getrusage isn't implemented by MS 00443 * also is not implemented on our HP superdome */ 00444 struct timeval { 00445 long tv_sec; /* seconds */ 00446 long tv_usec; /* microseconds */ 00447 }; 00448 #else 00449 #include <sys/time.h> 00450 #include <sys/resource.h> 00451 #endif 00452 00453 /* will be used to save initial time */ 00454 static struct timeval before; 00455 00456 /* cdClock stores time since arbitrary datum in clock_dat */ 00457 STATIC void cdClock(struct timeval *clock_dat) 00458 { 00459 DEBUG_ENTRY( "cdClock()" ); 00460 00461 /* >>chng 06 sep 2 rjrw: use long recurrence, fine grain UNIX clock * 00462 * -- maintain system dependences in a single routine */ 00463 #if defined(_MSC_VER) || defined(__HP_aCC) 00464 clock_t clock_val; 00465 /* >>chng 05 dec 21, from above to below due to negative exec times */ 00466 /*return( (double)(clock() - before) / (double)CLOCKS_PER_SEC );*/ 00467 clock_val = clock(); 00468 clock_dat->tv_sec = clock_val/CLOCKS_PER_SEC; 00469 clock_dat->tv_usec = 1000000*(clock_val-(clock_dat->tv_sec*CLOCKS_PER_SEC))/CLOCKS_PER_SEC; 00470 /*>>chng 06 oct 05, this produced very large number, time typically 50% too long 00471 clock_dat->tv_usec = 0;*/ 00472 #else 00473 struct rusage rusage; 00474 if(getrusage(RUSAGE_SELF,&rusage) != 0) 00475 { 00476 fprintf( ioQQQ, "DISASTER cdClock called getrusage with invalid arguments.\n" ); 00477 fprintf( ioQQQ, "Sorry.\n" ); 00478 cdEXIT(EXIT_FAILURE); 00479 } 00480 clock_dat->tv_sec = rusage.ru_utime.tv_sec; 00481 clock_dat->tv_usec = rusage.ru_utime.tv_usec; 00482 #endif 00483 00484 } 00485 /* cdSetExecTime is called by cdInit when everything is initialized, 00486 * so that every time cdExecTime is called the elapsed time is returned */ 00487 void cdSetExecTime(void) 00488 { 00489 cdClock(&before); 00490 lgCalled = true; 00491 } 00492 /* cdExecTime returns the elapsed time cpu time (sec) that has elapsed 00493 * since cdInit called cdSetExecTime.*/ 00494 double cdExecTime(void) 00495 { 00496 DEBUG_ENTRY( "cdExecTime()" ); 00497 00498 struct timeval clock_dat; 00499 /* check that we were properly initialized */ 00500 if( lgCalled) 00501 { 00502 cdClock(&clock_dat); 00503 /*fprintf(ioQQQ,"\n DEBUG sec %.2e usec %.2e\n", 00504 (double)(clock_dat.tv_sec-before.tv_sec), 00505 1e-6*(double)(clock_dat.tv_usec-before.tv_usec));*/ 00506 return (double)(clock_dat.tv_sec-before.tv_sec)+1e-6*(double)(clock_dat.tv_usec-before.tv_usec); 00507 } 00508 else 00509 { 00510 /* this is a big problem, we were asked for the elapsed time but 00511 * the timer was not initialized by calling SetExecTime */ 00512 fprintf( ioQQQ, "DISASTER cdExecTime was called before SetExecTime, impossible.\n" ); 00513 fprintf( ioQQQ, "Sorry.\n" ); 00514 cdEXIT(EXIT_FAILURE); 00515 } 00516 } 00517 00518 /************************************************************************* 00519 * 00520 * cdPrintCommands prints all input commands into file 00521 * 00522 ************************************************************************/ 00523 00524 /* cdPrintCommands( FILE *) 00525 * prints all input commands into file */ 00526 void cdPrintCommands( FILE * ioOUT ) 00527 { 00528 long int i; 00529 fprintf( ioOUT, " Input commands follow:\n" ); 00530 fprintf( ioOUT, "c ======================\n" ); 00531 00532 for( i=0; i <= input.nSave; i++ ) 00533 { 00534 fprintf( ioOUT, "%s\n", input.chCardSav[i] ); 00535 } 00536 fprintf( ioOUT, "c ======================\n" ); 00537 } 00538 00539 00540 /************************************************************************* 00541 * 00542 * cdEms obtain the local emissivity for a line, for the last computed zone 00543 * 00544 ************************************************************************/ 00545 00546 /* \todo 2 This routine, cdLine, cdEmis_ip, and cdLine_ip should be consolidated somewhat. 00547 * in particular so that this routine has the same "closest line" reporting that cdLine has. */ 00548 long int cdEmis( 00549 /* return value will be index of line within stack, 00550 * negative of number of lines in the stack if the line could not be found*/ 00551 /* 4 char null terminated string label */ 00552 char *chLabel, 00553 /* line wavelength */ 00554 realnum wavelength, 00555 /* the vol emissivity of this line in last computed zone */ 00556 double *emiss ) 00557 { 00558 /* use following to store local image of character strings */ 00559 char chCARD[INPUT_LINE_LENGTH]; 00560 char chCaps[5]; 00561 long int j; 00562 realnum errorwave; 00563 00564 DEBUG_ENTRY( "cdEms()" ); 00565 00566 /* routine returns the emissivity in the desired line 00567 * only used internally by the code, to do punch lines structure */ 00568 00569 strcpy( chCARD, chLabel ); 00570 00571 /* make sure chLabel is all caps */ 00572 caps(chCARD);/* convert to caps */ 00573 00574 /* get the error associated with 4 significant figures */ 00575 errorwave = WavlenErrorGet( wavelength ); 00576 00577 for( j=0; j < LineSave.nsum; j++ ) 00578 { 00579 /* change chLabel to all caps to be like input chLineLabel */ 00580 cap4(chCaps , (char*)LineSv[j].chALab); 00581 00582 /* check wavelength and chLabel for a match */ 00583 /*if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength)<errorwave && 00584 strcmp(chCaps,chCARD) == 0 ) */ 00585 if( fabs(LineSv[j].wavelength-wavelength) < errorwave && strcmp(chCaps,chCARD) == 0 ) 00586 { 00587 /* match, so set emiss to emissivity in line */ 00588 *emiss = LineSv[j].emslin[LineSave.lgLineEmergent]; 00589 /* and announce success by returning line index within stack */ 00590 return j; 00591 } 00592 } 00593 00594 /* we fell through without finding the line - return false */ 00595 return -LineSave.nsum; 00596 } 00597 00598 /************************************************************************* 00599 * 00600 * cdEms_ip obtain the local emissivity for a line with known index 00601 * 00602 ************************************************************************/ 00603 00604 void cdEmis_ip( 00605 /* index of the line in the stack */ 00606 long int ipLine, 00607 /* the vol emissivity of this line in last computed zone */ 00608 double *emiss ) 00609 { 00610 DEBUG_ENTRY( "cdEmis_ip()" ); 00611 00612 /* returns the emissivity in a line - implements punch lines structure 00613 * this version uses previously stored line index to save speed */ 00614 ASSERT( ipLine >= 0 && ipLine < LineSave.nsum ); 00615 *emiss = LineSv[ipLine].emslin[LineSave.lgLineEmergent]; 00616 return; 00617 } 00618 00619 /************************************************************************* 00620 * 00621 * cdColm get the column density for a constituent - 0 return if ok, 1 if problems 00622 * 00623 ************************************************************************/ 00624 00625 int cdColm( 00626 /* return value is zero if all ok, 1 if errors happened */ 00627 /* 4-char + eol string that is first 00628 * 4 char of element name as spelled by Cloudy, upper or lower case */ 00629 const char *chLabel, 00630 00631 /* integer stage of ionization, 1 for atom, 2 for A+, etc, 00632 * 0 is special flag for CO, H2, OH, or excited state */ 00633 long int ion, 00634 00635 /* the column density derived by the code [cm-2] */ 00636 double *theocl ) 00637 { 00638 long int nelem; 00639 /* use following to store local image of character strings */ 00640 char chLABEL_CAPS[20]; 00641 00642 DEBUG_ENTRY( "cdColm()" ); 00643 00644 /* check that chLabel[4] is null - supposed to be 4 char + end */ 00645 if( chLabel[4]!=0 ) 00646 { 00647 fprintf( ioQQQ, " cdColm called with insane ion label, =%s, must be 4 character + end of string.\n", 00648 chLabel ); 00649 return(1); 00650 } 00651 00652 strcpy( chLABEL_CAPS, chLabel ); 00653 /* convert element label to all caps */ 00654 caps(chLABEL_CAPS); 00655 00656 /* zero ionization stage has special meaning. The quantities recognized are 00657 * the molecules, "H2 ", "OH ", "CO ", etc 00658 * "CII*" excited state C+ */ 00659 if( ion < 0 ) 00660 { 00661 fprintf( ioQQQ, " cdColm called with insane ion, =%li\n", 00662 ion ); 00663 return(1); 00664 } 00665 00666 else if( ion == 0 ) 00667 { 00668 /* this case molecular column density */ 00669 /* want the molecular hydrogen H2 column density */ 00670 if( strcmp( chLABEL_CAPS , "H2 " )==0 ) 00671 { 00672 *theocl = colden.colden[ipCOL_H2g] + colden.colden[ipCOL_H2s]; 00673 } 00674 00675 /* H- column density */ 00676 else if( strcmp( chLABEL_CAPS , "H- " )==0 ) 00677 { 00678 *theocl = colden.colden[ipCOL_HMIN]; 00679 } 00680 00681 /* H2+ column density ipCOL_H2p is 4 */ 00682 else if( strcmp( chLABEL_CAPS , "H2+ " )==0 ) 00683 { 00684 *theocl = colden.colden[ipCOL_H2p]; 00685 } 00686 00687 /* H3+ column density */ 00688 else if( strcmp( chLABEL_CAPS , "H3+ " )==0 ) 00689 { 00690 *theocl = colden.colden[ipCOL_H3p]; 00691 } 00692 00693 /* H2g - ground H2 column density */ 00694 else if( strcmp( chLABEL_CAPS , "H2G " )==0 ) 00695 { 00696 *theocl = colden.colden[ipCOL_H2g]; 00697 } 00698 00699 /* H2* - excited H2 - column density */ 00700 else if( strcmp( chLABEL_CAPS , "H2* " )==0 ) 00701 { 00702 *theocl = colden.colden[ipCOL_H2s]; 00703 } 00704 00705 /* HeH+ column density */ 00706 else if( strcmp( chLABEL_CAPS , "HEH+" )==0 ) 00707 { 00708 *theocl = colden.colden[ipCOL_HeHp]; 00709 } 00710 00711 /* carbon monoxide column density */ 00712 else if( strcmp( chLABEL_CAPS , "CO " )==0 ) 00713 { 00714 *theocl = findspecies("CO")->hevcol; 00715 } 00716 00717 /* OH column density */ 00718 else if( strcmp( chLABEL_CAPS , "OH " )==0 ) 00719 { 00720 *theocl = findspecies("OH")->hevcol; 00721 } 00722 00723 /* H2O column density */ 00724 else if( strcmp( chLABEL_CAPS , "H2O " )==0 ) 00725 { 00726 *theocl = findspecies("H2O")->hevcol; 00727 } 00728 00729 /* O2 column density */ 00730 else if( strcmp( chLABEL_CAPS , "O2 " )==0 ) 00731 { 00732 *theocl = findspecies("O2")->hevcol; 00733 } 00734 00735 /* SiO column density */ 00736 else if( strcmp( chLABEL_CAPS , "SIO " )==0 ) 00737 { 00738 *theocl = findspecies("SiO")->hevcol; 00739 } 00740 00741 /* C2 column density */ 00742 else if( strcmp( chLABEL_CAPS , "C2 " )==0 ) 00743 { 00744 *theocl = findspecies("C2")->hevcol; 00745 } 00746 00747 /* C3 column density */ 00748 else if( strcmp( chLABEL_CAPS , "C3 " )==0 ) 00749 { 00750 *theocl = findspecies("C3")->hevcol; 00751 } 00752 00753 /* CN column density */ 00754 else if( strcmp( chLABEL_CAPS , "CN " )==0 ) 00755 { 00756 *theocl = findspecies("CN")->hevcol; 00757 } 00758 00759 /* CH column density */ 00760 else if( strcmp( chLABEL_CAPS , "CH " )==0 ) 00761 { 00762 *theocl = findspecies("CH")->hevcol; 00763 } 00764 00765 /* CH+ column density */ 00766 else if( strcmp( chLABEL_CAPS , "CH+ " )==0 ) 00767 { 00768 *theocl = findspecies("CH+")->hevcol; 00769 } 00770 00771 /* ===========================================================*/ 00772 /* end special case molecular column densities, start special cases 00773 * excited state column densities */ 00774 /* CII^* column density, population of J=3/2 upper level of split ground term */ 00775 else if( strcmp( chLABEL_CAPS , "CII*" )==0 ) 00776 { 00777 *theocl = colden.C2Colden[1]; 00778 } 00779 else if( strcmp( chLABEL_CAPS , "C11*" )==0 ) 00780 { 00781 *theocl = colden.C1Colden[0]; 00782 } 00783 else if( strcmp( chLABEL_CAPS , "C12*" )==0 ) 00784 { 00785 *theocl = colden.C1Colden[1]; 00786 } 00787 else if( strcmp( chLABEL_CAPS , "C13*" )==0 ) 00788 { 00789 *theocl = colden.C1Colden[2]; 00790 } 00791 else if( strcmp( chLABEL_CAPS , "O11*" )==0 ) 00792 { 00793 *theocl = colden.O1Colden[0]; 00794 } 00795 else if( strcmp( chLABEL_CAPS , "O12*" )==0 ) 00796 { 00797 *theocl = colden.O1Colden[1]; 00798 } 00799 else if( strcmp( chLABEL_CAPS , "O13*" )==0 ) 00800 { 00801 *theocl = colden.O1Colden[2]; 00802 } 00803 /* CIII excited states, upper level of 1909 */ 00804 else if( strcmp( chLABEL_CAPS , "C30*" )==0 ) 00805 { 00806 *theocl = colden.C3Colden[1]; 00807 } 00808 else if( strcmp( chLABEL_CAPS , "C31*" )==0 ) 00809 { 00810 *theocl = colden.C3Colden[2]; 00811 } 00812 else if( strcmp( chLABEL_CAPS , "C32*" )==0 ) 00813 { 00814 *theocl = colden.C3Colden[3]; 00815 } 00816 else if( strcmp( chLABEL_CAPS , "SI2*" )==0 ) 00817 { 00818 *theocl = colden.Si2Colden[1]; 00819 } 00820 else if( strcmp( chLABEL_CAPS , "HE1*" )==0 ) 00821 { 00822 *theocl = colden.He123S; 00823 } 00824 /* special option, "H2vJ" */ 00825 else if( strncmp(chLABEL_CAPS , "H2" , 2 ) == 0 ) 00826 { 00827 long int iVib = chLABEL_CAPS[2] - '0'; 00828 long int iRot = chLABEL_CAPS[3] - '0'; 00829 if( iVib<0 || iRot < 0 ) 00830 { 00831 fprintf( ioQQQ, " cdColm called with insane v,J for H2=\"%4.4s\" caps=\"%4.4s\"\n", 00832 chLabel , chLABEL_CAPS ); 00833 return( 1 ); 00834 } 00835 *theocl = cdH2_colden( iVib , iRot ); 00836 } 00837 00838 /* clueless as to what was meant - bail */ 00839 else 00840 { 00841 fprintf( ioQQQ, " cdColm called with unknown element chLabel, org=\"%4.4s \" 0 caps=\"%4.4s\" 0\n", 00842 chLabel , chLABEL_CAPS ); 00843 return(1); 00844 } 00845 } 00846 else 00847 { 00848 /* this case, ionization stage of some element */ 00849 /* find which element this is */ 00850 nelem = 0; 00851 while( nelem < LIMELM && 00852 strncmp(chLABEL_CAPS,elementnames.chElementNameShort[nelem],4) != 0 ) 00853 { 00854 ++nelem; 00855 } 00856 00857 /* this is true if we have one of the first 30 elements in the label, 00858 * nelem is on C scale */ 00859 if( nelem < LIMELM ) 00860 { 00861 00862 /* sanity check - does this ionization stage exist? 00863 * max2 is to pick up H2 as H 3 */ 00864 if( ion > MAX2(3,nelem + 2) ) 00865 { 00866 fprintf( ioQQQ, 00867 " cdColm asked to return ionization stage %ld for element %s but this is not physical.\n", 00868 ion, chLabel ); 00869 return(1); 00870 } 00871 00872 /* the column density, ion is on physics scale, but means are on C scale */ 00873 *theocl = mean.xIonMeans[0][nelem][ion-1]; 00874 /*>>chng 06 jan 23, div by factor of two 00875 * special case of H2 when being tricked as H 3 - this stores 2H_2 so that 00876 * the fraction of H in H0 and H+ is correct - need to remove this extra 00877 * factor of two here */ 00878 if( nelem==ipHYDROGEN && ion==3 ) 00879 *theocl /= 2.; 00880 } 00881 else 00882 { 00883 fprintf( ioQQQ, 00884 " cdColm did not understand this combination of ion %4ld and element %4.4s.\n", 00885 ion, chLabel ); 00886 return(1); 00887 } 00888 } 00889 return(0); 00890 } 00891 00892 00893 /************************************************************************* 00894 * 00895 *cdErrors produce summary of all warnings, cautions, etc, on arbitrary io unit 00896 * 00897 ************************************************************************/ 00898 00899 void cdErrors(FILE *ioOUT) 00900 { 00901 long int nc, 00902 nn, 00903 npe, 00904 ns, 00905 nte, 00906 nw , 00907 nIone, 00908 nEdene; 00909 bool lgAbort_loc; 00910 00911 DEBUG_ENTRY( "cdErrors()" ); 00912 00913 /* first check for number of warnings, cautions, etc */ 00914 cdNwcns(&lgAbort_loc,&nw,&nc,&nn,&ns,&nte,&npe, &nIone, &nEdene ); 00915 00916 /* only say something is one of these problems is nonzero */ 00917 if( nw || nc || nte || npe || nIone || nEdene || lgAbort_loc ) 00918 { 00919 /* say the title of the model */ 00920 fprintf( ioOUT, "%75.75s\n", input.chTitle ); 00921 00922 if( lgAbort_loc ) 00923 fprintf(ioOUT," Calculation ended with abort!\n"); 00924 00925 /* print warnings on the io unit */ 00926 if( nw != 0 ) 00927 { 00928 cdWarnings(ioOUT); 00929 } 00930 00931 /* print cautions on the io unit */ 00932 if( nc != 0 ) 00933 { 00934 cdCautions(ioOUT); 00935 } 00936 00937 if( nte != 0 ) 00938 { 00939 fprintf( ioOUT , "Te failures=%4ld\n", nte ); 00940 } 00941 00942 if( npe != 0 ) 00943 { 00944 fprintf( ioOUT , "Pressure failures=%4ld\n", npe ); 00945 } 00946 00947 if( nIone != 0 ) 00948 { 00949 fprintf( ioOUT , "Ionization failures=%4ld\n", nte ); 00950 } 00951 00952 if( nEdene != 0 ) 00953 { 00954 fprintf( ioOUT , "Electron density failures=%4ld\n", npe ); 00955 } 00956 } 00957 return; 00958 } 00959 00960 /************************************************************************* 00961 * 00962 *cdDepth_depth get depth structure from previous iteration 00963 * 00964 ************************************************************************/ 00965 void cdDepth_depth( double cdDepth[] ) 00966 { 00967 long int nz; 00968 00969 DEBUG_ENTRY( "cdDepth_depth()" ); 00970 00971 for( nz = 0; nz<nzone; ++nz ) 00972 { 00973 cdDepth[nz] = struc.depth[nz]; 00974 } 00975 return; 00976 } 00977 00978 /************************************************************************* 00979 * 00980 *cdPressure_depth routine to query results and return pressure of last iteration 00981 * 00982 ************************************************************************/ 00983 00984 /* 00985 * cdPressure_depth 00986 * This returns the pressure and its constituents for the last iteration. 00987 * space was allocated in the calling routine for the vectors - 00988 * before calling this, cdnZone should have been called to get the number of 00989 * zones, then space allocated for the arrays */ 00990 void cdPressure_depth( 00991 /* total pressure, all forms*/ 00992 double TotalPressure[], 00993 /* gas pressure */ 00994 double GasPressure[], 00995 /* line radiation pressure */ 00996 double RadiationPressure[]) 00997 { 00998 long int nz; 00999 01000 DEBUG_ENTRY( "cdPressure_depth()" ); 01001 01002 for( nz = 0; nz<nzone; ++nz ) 01003 { 01004 TotalPressure[nz] = struc.pressure[nz]; 01005 GasPressure[nz] = struc.GasPressure[nz]; 01006 RadiationPressure[nz] = struc.pres_radiation_lines_curr[nz]; 01007 } 01008 return; 01009 } 01010 01011 /************************************************************************* 01012 * 01013 *cdPressure_last routine to query results and return pressure of last zone 01014 * 01015 ************************************************************************/ 01016 01017 void cdPressure_last( 01018 double *PresTotal, /* total pressure, all forms, for the last computed zone*/ 01019 double *PresGas, /* gas pressure */ 01020 double *PresRad) /* line radiation pressure */ 01021 { 01022 01023 DEBUG_ENTRY( "cdPressure_last()" ); 01024 01025 *PresGas = pressure.PresGasCurr; 01026 *PresRad = pressure.pres_radiation_lines_curr; 01027 *PresTotal = pressure.PresTotlCurr; 01028 return; 01029 } 01030 01031 /************************************************************************* 01032 * 01033 *cdnZone gets number of zones 01034 * 01035 ************************************************************************/ 01036 01037 /* returns number of zones */ 01038 long int cdnZone( void ) 01039 { 01040 return nzone; 01041 } 01042 01043 /************************************************************************* 01044 * 01045 *cdTemp_last routine to query results and return temperature of last zone 01046 * 01047 ************************************************************************/ 01048 01049 01050 double cdTemp_last(void) 01051 { 01052 return phycon.te; 01053 } 01054 01055 /************************************************************************* 01056 * 01057 *cdIonFrac get ionization fractions for a constituent 01058 * 01059 ************************************************************************/ 01060 01061 int cdIonFrac( 01062 /* four char string, null terminated, giving the element name */ 01063 const char *chLabel, 01064 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 01065 * 0 says special case */ 01066 long int IonStage, 01067 /* will be fractional ionization */ 01068 double *fracin, 01069 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 01070 const char *chWeight , 01071 /* if true then weighting also has electron density, if false then only volume or radius */ 01072 bool lgDensity ) 01073 /* return value is 0 if element was found, non-zero if failed */ 01074 { 01075 01076 bool lgVol; 01077 long int ip, 01078 ion, /* used as index within aaa vector*/ 01079 nelem; 01080 realnum aaa[LIMELM + 1]; 01081 /* use following to store local image of character strings */ 01082 char chCARD[INPUT_LINE_LENGTH]; 01083 01084 DEBUG_ENTRY( "cdIonFrac()" ); 01085 01086 strcpy( chCARD, chWeight ); 01087 /* make sure chWeight is all caps */ 01088 caps(chCARD);/* convert to caps */ 01089 01090 /*caps(chWeight);*/ 01091 01092 if( strcmp(chCARD,"RADIUS") == 0 ) 01093 { 01094 lgVol = false; 01095 } 01096 01097 else if( strcmp(chCARD,"VOLUME") == 0 ) 01098 { 01099 lgVol = true; 01100 } 01101 01102 else 01103 { 01104 fprintf( ioQQQ, " cdIonFrac: chWeight=%6.6s makes no sense to me, valid options are RADIUS and VOLUME\n", 01105 chWeight ); 01106 *fracin = 0.; 01107 return(1); 01108 } 01109 01110 /* first ensure that chLabel is all caps */ 01111 strcpy( chCARD, chLabel ); 01112 /* make sure chLabel is all caps */ 01113 caps(chCARD);/* convert to caps */ 01114 01115 if( IonStage==0 ) 01116 { 01117 /* special case */ 01118 if( strcmp(chCARD,"H2 " ) == 0 ) 01119 { 01120 /* this will be request for H2, on c scale, hydro is 0 */ 01121 nelem = 0; 01122 IonStage = 3; 01123 } 01124 else 01125 { 01126 fprintf( ioQQQ, " cdIonFrac: ion stage of zero and element %s makes no sense to me\n", 01127 chCARD ); 01128 *fracin = 0.; 01129 return(1); 01130 } 01131 } 01132 01133 else 01134 { 01135 /* find which element this is, nelem is on physical scale, H is 1 */ 01136 nelem = 0; 01137 while( nelem < LIMELM && 01138 strcmp(chCARD,elementnames.chElementNameShort[nelem]) != 0 ) 01139 { 01140 ++nelem; 01141 } 01142 /* after this loop nelem is on c scale, H is 0 */ 01143 } 01144 01145 /* if element not recognized and above loop falls through, nelem is LIMELM 01146 * nelem counter is on physical scale, H = 1 Zn = 30 */ 01147 if( nelem >= LIMELM ) 01148 { 01149 fprintf( ioQQQ, " cdIonFrac called with unknown element chLabel, =%4.4s\n", 01150 chLabel ); 01151 return(1); 01152 } 01153 01154 /* sanity check - does this ionization stage exist? 01155 * both IonStage and nelem are on physical scales, H atoms are 1 1 */ 01156 01157 /* IonStage came in on physics scale, 01158 * ion will be used as pointer within the aaa array that contains average values, 01159 * convert to C scale */ 01160 ion = IonStage - 1; 01161 01162 /*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/ 01163 if( (ion > nelem+1 || ion < 0 ) && !(nelem==ipHYDROGEN&&ion==2)) 01164 { 01165 fprintf( ioQQQ, " cdIonFrac asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 01166 IonStage, chLabel ); 01167 *fracin = -1.; 01168 return(1); 01169 } 01170 01171 /* get either volume or radius average, aaa is filled in from 0 */ 01172 if( lgVol ) 01173 { 01174 /* aaa is dim'd LIMELM+1 so largest argument is LIMELM 01175 * 'i' means ionization, not temperature 01176 * nelem is on the physical scale */ 01177 /* last argument says whether to include electron density */ 01178 /* MeanIonVolume uses c scale for nelem */ 01179 MeanIonVolume('i', nelem,&ip,aaa,lgDensity); 01180 *fracin = pow((realnum)10.f,aaa[ion]); 01181 } 01182 else 01183 { 01184 /* last argument says whether to include electron density */ 01185 /* MeanIonRadius uses c scale for nelem */ 01186 MeanIonRadius('i', nelem,&ip,aaa,lgDensity); 01187 *fracin = pow((realnum)10.f,aaa[ion]); 01188 } 01189 01190 /* we succeeded - say so */ 01191 return(0); 01192 } 01193 01194 /************************************************************************* 01195 * 01196 * debugLine provides a debugging hook into the main line array 01197 * 01198 ************************************************************************/ 01199 01200 /* debugLine provides a debugging hook into the main line array 01201 * loops over whole array and finds every line that matches length, 01202 * the wavelength, the argument to the function 01203 * put breakpoint inside if test 01204 * return value is number of matches, also prints all matches*/ 01205 long debugLine( realnum wavelength ) 01206 { 01207 long j, kount; 01208 realnum errorwave; 01209 01210 kount = 0; 01211 01212 /* get the error associated with 4 significant figures */ 01213 errorwave = WavlenErrorGet( wavelength ); 01214 01215 for( j=0; j < LineSave.nsum; j++ ) 01216 { 01217 /* check wavelength and chLabel for a match */ 01218 /* if( fabs(LineSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave ) */ 01219 if( fabs(LineSv[j].wavelength-wavelength) < errorwave ) 01220 { 01221 printf("%s\n", LineSv[j].chALab); 01222 ++kount; 01223 } 01224 } 01225 printf(" hits = %li\n", kount ); 01226 return(kount); 01227 } 01228 01229 /************************************************************************* 01230 * 01231 *cdLineListPunch create a file with a list of all emission lines punched, 01232 *and their index within the emission line stack 01233 * 01234 ************************************************************************/ 01235 01236 /* returns total number of lines in the list */ 01237 long int cdLineListPunch( 01238 /* a file handle pointing to a file that is read for writing - 01239 * the calling routine must close it */ 01240 FILE* io ) 01241 { 01242 01243 long int j; 01244 01245 DEBUG_ENTRY( "cdLineListPunch()" ); 01246 01247 for( j=1; j < LineSave.nsum; j++ ) 01248 { 01249 01250 fprintf(io, "%li\t%s\t", 01251 j, 01252 LineSv[j].chALab); 01253 prt_wl( io , LineSv[j].wavelength ); 01254 fprintf(io, "\n"); 01255 } 01256 01257 01258 /* return total number of lines as debugging aid */ 01259 return LineSave.nsum; 01260 } 01261 01262 /************************************************************************* 01263 * 01264 *cdLine get the predicted line intensity, also index for line in stack 01265 * 01266 ************************************************************************/ 01267 01268 /* returns array index for line in array stack if we found the line, 01269 * return negative of total number of lines as debugging aid if line not found */ 01270 long int cdLine( 01271 const char *chLabel, 01272 /* wavelength of line in angstroms, not format printed by code */ 01273 realnum wavelength, 01274 /* linear intensity relative to normalization line*/ 01275 double *relint, 01276 /* log of luminosity or intensity of line */ 01277 double *absint ) 01278 { 01279 char chCaps[5], 01280 chFind[5]; 01281 long int ipobs, 01282 /* >>chng 05 dec 21, RP add option to print closest line when 01283 * we don't find the line we are after */ 01284 j, index_of_closest=LONG_MIN, 01285 index_of_closest_w_correct_label=-1; 01286 realnum errorwave, smallest_error=BIGFLOAT, 01287 smallest_error_w_correct_label=BIGFLOAT; 01288 char chLabelLoc[INPUT_LINE_LENGTH]; 01289 01290 DEBUG_ENTRY( "cdLine()" ); 01291 01292 /* this is zero when cdLine called with cdNoExec called too */ 01293 if( LineSave.nsum == 0 ) 01294 { 01295 *relint = 0.; 01296 *absint = 0.; 01297 return 0; 01298 } 01299 ASSERT( LineSave.ipNormWavL >= 0); 01300 ASSERT( LineSave.nsum > 0); 01301 01302 /* check that chLabel[4] is null - supposed to be 4 char + end */ 01303 if( chLabel[4]!=0 ) 01304 { 01305 fprintf( ioQQQ, " cdLine called with insane line label, =%s, must be 4 character + end of string.\n", 01306 chLabel ); 01307 return 1; 01308 } 01309 01310 /* change chLabel to all caps */ 01311 strcpy( chLabelLoc , chLabel ); 01312 /*cap4(chFind,chLabel);*/ 01313 cap4(chFind,chLabelLoc); 01314 01315 /* this replaces tabs with spaces. */ 01316 /* \todo 2 keep this in, do it elsewhere, just warn and bail? */ 01317 for( j=0; j<=3; j++ ) 01318 { 01319 if( chFind[j] == '\t' ) 01320 { 01321 chFind[j] = ' '; 01322 } 01323 } 01324 01325 /* get the error associated with 4 significant figures */ 01326 errorwave = WavlenErrorGet( wavelength ); 01327 01328 { 01329 /*@-redef@*/ 01330 enum{DEBUG_LOC=false}; 01331 /*@+redef@*/ 01332 if( DEBUG_LOC && fabs(wavelength-1000.)<0.01 ) 01333 { 01334 fprintf(ioQQQ,"cdDrive wl %.4e error %.3e\n", 01335 wavelength, errorwave ); 01336 } 01337 } 01338 01339 /* now go through entire line stack, do not do 0, which is unity integration */ 01340 for( j=1; j < LineSave.nsum; j++ ) 01341 { 01342 /* >>chng 05 dec 21, find closest line to requested wavelength to 01343 * report when we don't get exact match */ 01344 realnum current_error; 01345 current_error = (realnum)fabs(LineSv[j].wavelength-wavelength); 01346 01347 /* change chLabel to all caps to be like input chALab */ 01348 cap4(chCaps , (char*)LineSv[j].chALab); 01349 01350 if( current_error < smallest_error ) 01351 { 01352 index_of_closest = j; 01353 smallest_error = current_error; 01354 } 01355 01356 if( current_error < smallest_error_w_correct_label && 01357 (strcmp(chCaps,chFind) == 0) ) 01358 { 01359 index_of_closest_w_correct_label = j; 01360 smallest_error_w_correct_label = current_error; 01361 } 01362 01363 /* check wavelength and chLabel for a match */ 01364 /* DELTA since often want wavelength of zero */ 01365 /*if( fabs(LineSv[j].wavelength-wavelength)/MAX2(DELTA,wavelength) < errorwave )*/ 01366 if( current_error < errorwave ) 01367 { 01368 01369 /* change chLabel to all caps to be like input chALab 01370 cap4(chCaps , (char*)LineSv[j].chALab);*/ 01371 01372 /* now see if labels agree */ 01373 if( strcmp(chCaps,chFind) == 0 ) 01374 { 01375 /* match, so set pointer */ 01376 ipobs = j; 01377 01378 /* does the normalization line have a positive intensity*/ 01379 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 01380 { 01381 *relint = LineSv[ipobs].sumlin[LineSave.lgLineEmergent]/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]* 01382 LineSave.ScaleNormLine; 01383 } 01384 else 01385 { 01386 *relint = 0.; 01387 } 01388 01389 /* return log of current line intensity if it is positive */ 01390 if( LineSv[ipobs].sumlin[LineSave.lgLineEmergent] > 0. ) 01391 { 01392 *absint = log10(LineSv[ipobs].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten; 01393 } 01394 else 01395 { 01396 /* line intensity is actually zero, return small number */ 01397 *absint = -37.; 01398 } 01399 /* we found the line, return pointer to its location */ 01400 return ipobs; 01401 } 01402 } 01403 } 01404 01405 /* >>chng 05 dec 21, report closest line if we did not find exact match, note that 01406 * exact match returns above, where we will return negative number of lines in stack */ 01407 fprintf( ioQQQ," NOTE cdLine did not find line with label %4s and wavelength ", chFind ); 01408 prt_wl(ioQQQ,wavelength); 01409 if( index_of_closest >= 0 ) 01410 { 01411 fprintf( ioQQQ,".\n The closest line (any label) was \t%4s\t", LineSv[index_of_closest].chALab ); 01412 prt_wl(ioQQQ,LineSv[index_of_closest].wavelength ); 01413 /* fprintf( ioQQQ," %f", LineSv[index_of_closest].wavelength ); */ 01414 fprintf( ioQQQ,"\n The closest with correct label was\t%4s\t", chFind ); 01415 prt_wl(ioQQQ,LineSv[index_of_closest_w_correct_label].wavelength ); 01416 /*fprintf( ioQQQ," %f", LineSv[index_of_closest_w_correct_label].wavelength ); */ 01417 fprintf( ioQQQ,"\n" ); 01418 } 01419 else 01420 { 01421 fprintf( ioQQQ,".\n PROBLEM No close line was found\n" ); 01422 TotalInsanity(); 01423 } 01424 01425 *absint = 0.; 01426 *relint = 0.; 01427 01428 /* if we fell down to here we did not find the line */ 01429 /* >>chng 00 sep 02, return negative of total number of lines as debugging aid */ 01430 return -LineSave.nsum; 01431 } 01432 01433 /*cdLine_ip get the predicted line intensity, using index for line in stack */ 01434 void cdLine_ip(long int ipLine, 01435 /* linear intensity relative to normalization line*/ 01436 double *relint, 01437 /* log of luminosity or intensity of line */ 01438 double *absint ) 01439 { 01440 01441 DEBUG_ENTRY( "cdLine_ip()" ); 01442 01443 /* this is zero when cdLine called with cdNoExec called too */ 01444 if( LineSave.nsum == 0 ) 01445 { 01446 *relint = 0.; 01447 *absint = 0.; 01448 return; 01449 } 01450 ASSERT( LineSave.ipNormWavL >= 0); 01451 ASSERT( LineSave.nsum > 0); 01452 01453 /* does the normalization line have a positive intensity*/ 01454 if( LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent] > 0. ) 01455 { 01456 *relint = LineSv[ipLine].sumlin[LineSave.lgLineEmergent]/LineSv[LineSave.ipNormWavL].sumlin[LineSave.lgLineEmergent]* 01457 LineSave.ScaleNormLine; 01458 } 01459 else 01460 { 01461 *relint = 0.; 01462 } 01463 01464 /* return log of current line intensity if it is positive */ 01465 if( LineSv[ipLine].sumlin[LineSave.lgLineEmergent] > 0. ) 01466 { 01467 *absint = log10(LineSv[ipLine].sumlin[LineSave.lgLineEmergent]) + radius.Conv2PrtInten; 01468 } 01469 else 01470 { 01471 /* line intensity is actually zero, return small number */ 01472 *absint = -37.; 01473 } 01474 01475 return; 01476 } 01477 01478 /************************************************************************* 01479 * 01480 *cdDLine get the predicted emergent line intensity, also index for line in stack 01481 * 01482 ************************************************************************/ 01483 01484 /* returns array index for line in array stack if we found the line, 01485 * or false==0 if we did not find the line */ 01486 long int cdDLine(char *chLabel, 01487 /* wavelength of line as printed by code*/ 01488 realnum wavelength, 01489 /* linear intensity relative to normalization line*/ 01490 double *relint, 01491 /* log of luminosity or intensity of line */ 01492 double *absint ) 01493 { 01494 char chCaps[5], 01495 chFind[5]; 01496 long int ipobs, 01497 j; 01498 realnum errorwave; 01499 01500 DEBUG_ENTRY( "cdDLine()" ); 01501 01502 /* this is zero when cdDLine called with cdNoExec called too */ 01503 if( LineSave.nsum == 0 ) 01504 { 01505 *relint = 0.; 01506 *absint = 0.; 01507 return 0; 01508 } 01509 ASSERT( LineSave.ipNormWavL >= 0); 01510 ASSERT( LineSave.nsum > 0); 01511 01512 /* change chLabel to all caps */ 01513 cap4(chFind,chLabel); 01514 01515 /* get the error associated with 4 significant figures */ 01516 errorwave = WavlenErrorGet( wavelength ); 01517 01518 /* now go through entire line stack, do not do 0, which is H-beta and 01519 * in stack further down - this is to stop query for H-beta from returning 01520 * 0, the flag for line not found */ 01521 /* >>chng 06 mar 09, move to common line array */ 01522 for( j=1; j < LineSave.nsum; j++ ) 01523 { 01524 01525 /* check wavelength and chLabel for a match */ 01526 /* DELTA since often want wavelength of zero */ 01527 /* if( fabs(LineDSv[j].wavelength- wavelength)/MAX2(DELTA,wavelength) < errorwave ) */ 01528 /* >>chng 06 mar 09, first wavelength had been LineSv, should have been LineDSv, 01529 * so this routine was catching the wrong line */ 01530 if( fabs(LineSv[j].wavelength-wavelength) < errorwave ) 01531 { 01532 01533 /* change chLabel to all caps to be like input chALab */ 01534 cap4(chCaps , (char*)LineSv[j].chALab); 01535 01536 /* now see if labels agree */ 01537 if( strcmp(chCaps,chFind) == 0 ) 01538 { 01539 /* match, so set array index */ 01540 ipobs = j; 01541 01542 /* does the normalization line have a positive intensity*/ 01543 if( LineSv[LineSave.ipNormWavL].sumlin[1] > 0. ) 01544 { 01545 *relint = LineSv[ipobs].sumlin[1]/LineSv[LineSave.ipNormWavL].sumlin[1]* 01546 LineSave.ScaleNormLine; 01547 } 01548 else 01549 { 01550 *relint = 0.; 01551 } 01552 01553 /* return log of current line intensity if it is positive */ 01554 if( LineSv[ipobs].sumlin[1] > 0. ) 01555 { 01556 *absint = log10(LineSv[ipobs].sumlin[1]) + radius.Conv2PrtInten; 01557 } 01558 else 01559 { 01560 /* line intensity is actually zero, return small number */ 01561 *absint = -37.; 01562 } 01563 /* we found the line, return pointer to its location */ 01564 return ipobs; 01565 } 01566 } 01567 } 01568 01569 *absint = 0.; 01570 *relint = 0.; 01571 01572 /* if we fell down to here we did not find the line */ 01573 /* >>chng 00 sep 02, return negative of total number of lines as debugging aid */ 01574 return -LineSave.nsum; 01575 } 01576 01577 /************************************************************************* 01578 * 01579 *cdNwcns get the number of cautions and warnings, to tell if calculation is ok 01580 * 01581 ************************************************************************/ 01582 01583 void cdNwcns( 01584 /* abort status, this better be false */ 01585 bool *lgAbort_ret , 01586 /* the number of warnings, cautions, notes, and surprises */ 01587 long int *NumberWarnings, 01588 long int *NumberCautions, 01589 long int *NumberNotes, 01590 long int *NumberSurprises, 01591 /* the number of temperature convergence failures */ 01592 long int *NumberTempFailures, 01593 /* the number of pressure convergence failures */ 01594 long int *NumberPresFailures, 01595 /* the number of ionization convergence failures */ 01596 long int *NumberIonFailures, 01597 /* the number of electron density convergence failures */ 01598 long int *NumberNeFailures ) 01599 { 01600 01601 DEBUG_ENTRY( "cdNwcns()" ); 01602 01603 /* this would be set true if code ended with abort - very very bad */ 01604 *lgAbort_ret = lgAbort; 01605 /* this sub is called after comment, to find the number of various comments */ 01606 *NumberWarnings = warnings.nwarn; 01607 *NumberCautions = warnings.ncaun; 01608 *NumberNotes = warnings.nnote; 01609 *NumberSurprises = warnings.nbang; 01610 01611 /* these are counters that were incremented during convergence failures */ 01612 *NumberTempFailures = conv.nTeFail; 01613 *NumberPresFailures = conv.nPreFail; 01614 *NumberIonFailures = conv.nIonFail; 01615 *NumberNeFailures = conv.nNeFail; 01616 return; 01617 } 01618 01619 /************************************************************************* 01620 * 01621 * cdOutp redirect output to arbitrary opened file 01622 * 01623 ************************************************************************/ 01624 01625 void cdOutp( FILE *ioOut) 01626 { 01627 01628 DEBUG_ENTRY( "cdOutp()" ); 01629 01630 /* ioQQQ is pointer to output fiile */ 01631 ioQQQ = ioOut; 01632 return; 01633 } 01634 01635 /************************************************************************* 01636 * 01637 * cdInp redirect line input from arbitrary opened file 01638 * 01639 ************************************************************************/ 01640 01641 void cdInp( FILE *ioInp) 01642 { 01643 01644 DEBUG_ENTRY( "cdInp()" ); 01645 01646 /* ioQQQ is pointer to output file */ 01647 ioStdin = ioInp; 01648 return; 01649 } 01650 01651 /************************************************************************* 01652 * 01653 *cdTalk tells the code whether to print results or be silent 01654 * 01655 ************************************************************************/ 01656 01657 void cdTalk(bool lgTOn) 01658 { 01659 01660 DEBUG_ENTRY( "cdTalk()" ); 01661 01662 called.lgTalk = lgTOn; 01663 01664 if( lgTOn ) 01665 { 01666 /* means talk has not been forced off */ 01667 called.lgTalkForcedOff = false; 01668 } 01669 else 01670 { 01671 /* means talk is forced off */ 01672 called.lgTalkForcedOff = true; 01673 } 01674 return; 01675 } 01676 01677 /*cdB21cm - returns B as measured by 21 cm */ 01678 double cdB21cm( void ) 01679 { 01680 double ret; 01681 01682 DEBUG_ENTRY( "cdB21cm()" ); 01683 01684 if( mean.B_HarMeanTempRadius[1] > SMALLFLOAT ) 01685 { 01686 ret = mean.B_HarMeanTempRadius[0]/mean.B_HarMeanTempRadius[1]; 01687 } 01688 else 01689 { 01690 ret = 0.; 01691 } 01692 return(ret); 01693 } 01694 01695 /************************************************************************* 01696 * 01697 * cdTemp get mean electron temperature for any element 01698 * 01699 ************************************************************************/ 01700 01701 /* This routine finds the mean electron temperature for any ionization stage 01702 * It returns 0 if it could find the species, 1 if it could not find the species. 01703 * The first argument is a null terminated 4 char string that gives the element 01704 * name as spelled by Cloudy. 01705 * The second argument is ion stage, 1 for atom, 2 for first ion, etc 01706 * This third argument will be returned as result, 01707 * Last parameter is either "VOLUME" or "RADIUS" to give weighting 01708 * 01709 * if the ion stage is zero then the element label will have a special meaning. 01710 * The string "21CM" is will return the 21 cm temperature. 01711 * The string "H2 " will return the temperature weighted wrt the H2 density */ 01715 /* return value is o if things are ok and element was found, 01716 * non-zero if element not found or there are problems */ 01717 int cdTemp( 01718 /* four char string, null terminated, giving the element name */ 01719 const char *chLabel, 01720 /* IonStage is ionization stage, 1 for atom, up to N+1 where N is atomic number, 01721 * 0 means that chLabel is a special case */ 01722 long int IonStage, 01723 /* will be temperature */ 01724 double *TeMean, 01725 /* how to weight the average, must be "VOLUME" or "RADIUS" */ 01726 const char *chWeight ) 01727 { 01728 01729 bool lgVol; 01730 long int ip, 01731 ion, /* used as pointer within aaa vector*/ 01732 nelem; 01733 realnum aaa[LIMELM + 1]; 01734 /* use following to store local image of character strings */ 01735 char chWGHT[INPUT_LINE_LENGTH] , chELEM[INPUT_LINE_LENGTH]; 01736 01737 DEBUG_ENTRY( "cdTemp()" ); 01738 01739 /* make sure chWeight is all caps */ 01740 strcpy( chWGHT, chWeight ); 01741 caps(chWGHT);/* convert to caps */ 01742 01743 /* ensure that chLabel is all caps */ 01744 strcpy( chELEM, chLabel ); 01745 caps(chELEM);/* convert to caps */ 01746 01747 /* now see if volume or radius weighting */ 01748 if( strcmp(chWGHT,"RADIUS") == 0 ) 01749 { 01750 lgVol = false; 01751 } 01752 else if( strcmp(chWGHT,"VOLUME") == 0 ) 01753 { 01754 lgVol = true; 01755 } 01756 else 01757 { 01758 fprintf( ioQQQ, " cdTemp: chWeight=%6.6s makes no sense to me, the options are RADIUS and VOLUME.\n", 01759 chWeight ); 01760 *TeMean = 0.; 01761 return(1); 01762 } 01763 01764 if( IonStage == 0 ) 01765 { 01766 /* return atomic hydrogen weighted harmonic mean of gas kinetic temperature */ 01767 if( strcmp(chELEM,"21CM") == 0 ) 01768 { 01769 if( lgVol ) 01770 { 01771 /* this is mean of inverse temperature weighted over volume */ 01772 if( mean.HarMeanTempVolume[1] > SMALLFLOAT ) 01773 { 01774 *TeMean = mean.HarMeanTempVolume[0]/mean.HarMeanTempVolume[1]; 01775 } 01776 else 01777 { 01778 *TeMean = 0.; 01779 } 01780 } 01781 else 01782 { 01783 /* this is mean of inverse temperature weighted over radius */ 01784 if( mean.HarMeanTempRadius[1] > SMALLFLOAT ) 01785 { 01786 *TeMean = mean.HarMeanTempRadius[0]/mean.HarMeanTempRadius[1]; 01787 } 01788 else 01789 { 01790 *TeMean = 0.; 01791 } 01792 } 01793 } 01794 /* return atomic hydrogen weighted harmonic mean 21 cm spin temperature, 01795 * using actual level populations with 1s of H0 */ 01796 else if( strcmp(chELEM,"SPIN") == 0 ) 01797 { 01798 *TeMean = mean.H_21cm_spin_mean_radius[0] / 01799 SDIV(mean.H_21cm_spin_mean_radius[1]); 01800 } 01801 /* return temperature deduced from ratio of 21 cm and Lya optical depths */ 01802 else if( strcmp(chELEM,"OPTI") == 0 ) 01803 { 01804 /* this is the temperature derived from Lya - 21 cm optical depths */ 01805 *TeMean = 01806 3.84e-7 * Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon / 01807 SDIV( HFLines[0].Emis->TauCon ); 01808 } 01809 /* special case, mean wrt H_2 */ 01810 else if( strcmp(chELEM,"H2 ") == 0 ) 01811 { 01812 if( lgVol ) 01813 { 01814 /* mean temp with H2 over volume */ 01815 if( mean.H2MeanTempVolume[1] > SMALLFLOAT ) 01816 { 01817 *TeMean = mean.H2MeanTempVolume[0] / mean.H2MeanTempVolume[1]; 01818 } 01819 else 01820 { 01821 *TeMean = 0.; 01822 } 01823 } 01824 else 01825 { 01826 /* mean temp with H2 over radius */ 01827 if( mean.H2MeanTempRadius[1] > SMALLFLOAT ) 01828 { 01829 *TeMean = mean.H2MeanTempRadius[0] / mean.H2MeanTempRadius[1]; 01830 } 01831 else 01832 { 01833 *TeMean = 0.; 01834 } 01835 } 01836 } 01837 /* four spaces mean to return simple mean over rad or vol */ 01838 else if( strcmp(chELEM," ") == 0 ) 01839 { 01840 if( lgVol ) 01841 { 01842 /* this is temperature weighted over volume */ 01843 if( mean.TempMeanVolume[1] > SMALLFLOAT ) 01844 { 01845 *TeMean = mean.TempMeanVolume[0]/mean.TempMeanVolume[1]; 01846 } 01847 else 01848 { 01849 *TeMean = 0.; 01850 } 01851 } 01852 else 01853 { 01854 /* this is mean of inverse temperature weighted over radius */ 01855 if( mean.TempMeanRadius[1] > SMALLFLOAT ) 01856 { 01857 *TeMean = mean.TempMeanRadius[0]/mean.TempMeanRadius[1]; 01858 } 01859 else 01860 { 01861 *TeMean = 0.; 01862 } 01863 } 01864 } 01865 else 01866 { 01867 fprintf( ioQQQ, " cdTemp called with ion=0 and unknown quantity, =%4.4s\n", 01868 chLabel ); 01869 fprintf( ioQQQ, " I know about 21CM and H2__ \n"); 01870 return(1); 01871 } 01872 01873 /* say things are ok */ 01874 return(0); 01875 } 01876 01877 /* find which element this is */ 01878 nelem = 0; 01879 while( nelem < LIMELM && 01880 strcmp(chELEM,elementnames.chElementNameShort[nelem]) != 0 ) 01881 { 01882 ++nelem; 01883 } 01884 /* after this loop nelem is atomic number of element, H is 1 */ 01885 01886 /* if element not recognized and above loop falls through, nelem is LIMELM+1 01887 * nelem counter is on physical scale, H = 1 Zn = 30 */ 01888 if( nelem >= LIMELM ) 01889 { 01890 fprintf( ioQQQ, " cdTemp called with unknown element chLabel, =%4.4s\n", 01891 chLabel ); 01892 return(1); 01893 } 01894 01895 /* sanity check - does this ionization stage exist? 01896 * both IonStage on physical scale, nelem on c */ 01897 01898 /* ion will be used as pointer within the aaa array that contains average values, 01899 * done this way to prevent lint from false problem in access to aaa array */ 01900 ion = IonStage - 1; 01901 01902 /*if( IonStage > nelem + 1 || IonStage < 1 || IonStage > LIMELM+1 )*/ 01903 if( ion > nelem+1 || ion < 0 || ion > LIMELM ) 01904 { 01905 fprintf( ioQQQ, " cdTemp asked to return ionization stage%4ld for element %4.4s but this is not physical.\n", 01906 IonStage, chLabel ); 01907 return(1); 01908 } 01909 01910 /* get either volume or radius average, aaa is filled in from 0 */ 01911 if( lgVol ) 01912 { 01913 /* aaa is dim'd LIMELM+1 so largest arg is LIMELM */ 01914 /* MeanIonVolume uses C scale for nelem */ 01915 MeanIonVolume('t', nelem,&ip,aaa,false); 01916 *TeMean = pow((realnum)10.f,aaa[ion]); 01917 } 01918 else 01919 { 01920 MeanIonRadius('t', nelem,&ip,aaa,false); 01921 *TeMean = pow((realnum)10.f,aaa[ion]); 01922 } 01923 return(0); 01924 } 01925 01926 /************************************************************************* 01927 * 01928 * cdRead routine to read in command lines 01929 * called by user when cloudy used as subroutine 01930 * called by maincl when used as a routine 01931 * 01932 ************************************************************************/ 01933 01934 /* returns the number of lines available in command stack 01935 * this is limit to how many more commands can be read */ 01936 int cdRead( 01937 /* the string containing the command */ 01938 const char *chInputLine ) 01939 { 01940 char *chEOL , /* will be used to search for end of line symbols */ 01941 chCARD[INPUT_LINE_LENGTH], 01942 chLocal[INPUT_LINE_LENGTH]; 01943 01944 DEBUG_ENTRY( "cdRead()" ); 01945 01946 /* this is set false when code loaded, set true when cdInit called, 01947 * this is check that cdInit was called first */ 01948 if( !lgcdInitCalled ) 01949 { 01950 printf(" cdInit was not called first - this must be the first call.\n"); 01951 cdEXIT(EXIT_FAILURE); 01952 } 01953 01954 /* totally ignore any card starting with a #, *, space, //, or % 01955 * but want to include special "c " type of comment 01956 * >>chng 06 sep 04 use routine to check for comments */ 01957 if( ( lgInputComment( chInputLine ) || 01958 /* these two are end-of-input-stream sentinels */ 01959 chInputLine[0] == '\n' || chInputLine[0] == ' ' ) && 01960 /* option to allow "C" lines through */ 01961 ! ( chInputLine[0] == 'C' || chInputLine[0] == 'c' ) ) 01962 { 01963 /* return value is number of lines that can still be stuffed in */ 01964 return(NKRD - input.nSave); 01965 } 01966 01967 /*************************************************************************** 01968 * validate a location to store this line image, then store the version * 01969 * that has been truncated from special end of line characters * 01970 * stored image will have < INPUT_LINE_LENGTH valid characters * 01971 ***************************************************************************/ 01972 01973 /* this will now point to the next free slot in the line image save array 01974 * this is where we will stuff this line image */ 01975 ++input.nSave; 01976 01977 if( input.nSave >= NKRD ) 01978 { 01979 /* too many input commands were entered - bail */ 01980 fprintf( ioQQQ, 01981 " Too many line images entered to cdRead. The limit is %d\n", 01982 NKRD ); 01983 fprintf( ioQQQ, 01984 " The limit to the number of allowed input lines is %d. This limit was exceeded. Sorry.\n", 01985 NKRD ); 01986 fprintf( ioQQQ, 01987 " This limit is set by the variable NKRD which appears in input.h \n" ); 01988 fprintf( ioQQQ, 01989 " Increase it everywhere it appears.\n" ); 01990 cdEXIT(EXIT_FAILURE); 01991 } 01992 01993 strncpy( chLocal, chInputLine, INPUT_LINE_LENGTH ); 01994 // strncpy will pad chLocal with null bytes if chInputLine is shorter than 01995 // INPUT_LINE_LENGTH characters, so this indicates an overlong input string 01996 if( chLocal[INPUT_LINE_LENGTH-1] != '\0' ) 01997 { 01998 chLocal[INPUT_LINE_LENGTH-1] = '\0'; 01999 fprintf(ioQQQ," PROBLEM cdRead, while parsing the following input line:\n %s\n", 02000 chInputLine); 02001 fprintf(ioQQQ," found that the line is longer than %i characters, the longest possible line.\n", 02002 INPUT_LINE_LENGTH-1); 02003 fprintf(ioQQQ," Please make the command line no longer than this limit.\n"); 02004 } 02005 02006 /* now kill any part of line image after special end of line character, 02007 * this stops info user wants ignored from entering after here */ 02008 if( (chEOL = strchr(chLocal, '\n' ) ) != NULL ) 02009 { 02010 *chEOL = '\0'; 02011 } 02012 if( (chEOL = strchr(chLocal, '%' ) ) != NULL ) 02013 { 02014 *chEOL = '\0'; 02015 } 02016 /* >>chng 02 apr 10, add this char */ 02017 if( (chEOL = strchr(chLocal, '#' ) ) != NULL ) 02018 { 02019 *chEOL = '\0'; 02020 } 02021 if( (chEOL = strchr(chLocal, ';' ) ) != NULL ) 02022 { 02023 *chEOL = '\0'; 02024 } 02025 if( (chEOL = strstr(chLocal, "//" ) ) != NULL ) 02026 { 02027 *chEOL = '\0'; 02028 } 02029 02030 /* now do it again, since we now want to make sure that there is a trailing space 02031 * if the line is shorter than 80 char, test on null is to keep lint happy */ 02032 if( (chEOL = strchr( chLocal, '\0' )) == NULL ) 02033 TotalInsanity(); 02034 02035 /* pad with two spaces if short enough, 02036 * if not short enough for this to be done, then up to user to create correct input */ 02037 if( chEOL-chLocal < INPUT_LINE_LENGTH-2 ) 02038 { 02039 strcat( chLocal, " " ); 02040 } 02041 02042 /* save string in master array for later use in readr */ 02043 strcpy( input.chCardSav[input.nSave], chLocal ); 02044 02045 /* copy string over the chCARD, then convert this to caps to check for optimize*/ 02046 strcpy( chCARD, chLocal ); 02047 caps(chCARD);/* convert to caps */ 02048 02049 /* check whether this is a trace command, turn on printout if so */ 02050 if( strncmp(chCARD,"TRACE",5) == 0 ) 02051 { 02052 trace.lgTrace = true; 02053 } 02054 02055 /* print input lines if trace specified */ 02056 if( trace.lgTrace ) 02057 { 02058 fprintf( ioQQQ,"cdRead=%s=\n",input.chCardSav[input.nSave] ); 02059 } 02060 02061 /* now check whether VARY is specified */ 02062 if( nMatch("VARY",chCARD) ) 02063 { 02064 /* a optimize flag was on the line image */ 02065 optimize.lgVaryOn = true; 02066 } 02067 02068 /* now check whether line is "no vary" command */ 02069 if( strncmp(chCARD,"NO VARY",7) == 0 ) 02070 { 02071 optimize.lgNoVary = true; 02072 } 02073 02074 /* now check whether line is "grid" command */ 02075 if( strncmp(chCARD,"GRID",4) == 0 ) 02076 { 02077 grid.lgGrid = true; 02078 ++grid.nGridCommands; 02079 } 02080 02081 /* NO BUFFERING turn off buffered io for standard output, 02082 * used to get complete output when code crashes */ 02083 if( strncmp(chCARD,"NO BUFF",7) == 0 ) 02084 { 02085 /* if output has already been redirected (e.g. using cdOutp()) then 02086 * ignore this command, a warning will be printed in ParseDont() */ 02087 /* stdout may be a preprocessor macro, so lets be really careful here */ 02088 FILE *test = stdout; 02089 if( ioQQQ == test ) 02090 { 02091 fprintf( ioQQQ, " cdRead found NO BUFFERING command, redirecting output to stderr now.\n" ); 02092 /* make sure output is not lost */ 02093 fflush( ioQQQ ); 02094 /* stderr is always unbuffered */ 02095 ioQQQ = stderr; 02096 /* will be used to generate comment at end */ 02097 input.lgSetNoBuffering = false; 02098 } 02099 } 02100 02101 /* >>chng 06 may 20, use grid command as substitute for optimize command */ 02102 if( strncmp(chCARD,"OPTI",4) == 0 || strncmp(chCARD,"GRID",4) == 0 ) 02103 { 02104 /* optimize command read in */ 02105 optimize.lgOptimr = true; 02106 } 02107 return( NKRD - input.nSave ); 02108 } 02109 02110 /* wrapper to close all punch files */ 02111 void cdClosePunchFiles( void ) 02112 { 02113 02114 DEBUG_ENTRY( "cdClosePunchFiles()" ); 02115 02116 ClosePunchFiles( true ); 02117 return; 02118 } 02119 02120 #if 0 02121 /*cdTemp21cmLya this is the temperature derived from Lya - 21 cm optical depths */ 02122 double cdTemp21cmLya( void ) 02123 { 02124 double t; 02125 02126 DEBUG_ENTRY( "cdTemp21cmLya()" ); 02127 02128 /*fprintf(ioQQQ,"\n bug debug 21cm %.3e %.3e %.3e\n", 02129 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon , 02130 HFLines[0].TauCon , 02131 3.84e-7 * Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon / 02132 SDIV( HFLines[0].TauCon ) );*/ 02133 02134 /* this is the temperature derived from Lya - 21 cm optical depths */ 02135 t = 02136 3.84e-7 * Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon / 02137 SDIV( HFLines[0].TauCon ); 02138 return t; 02139 } 02140 #endif 02141