00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "cddefines.h"
00041 #include "trace.h"
00042 #include "conv.h"
00043 #include "init.h"
00044 #include "lines.h"
00045 #include "pressure.h"
00046 #include "prt.h"
00047 #include "colden.h"
00048 #include "dense.h"
00049 #include "radius.h"
00050 #include "struc.h"
00051 #include "mole.h"
00052 #include "elementnames.h"
00053 #include "mean.h"
00054 #include "phycon.h"
00055 #include "called.h"
00056 #include "parse.h"
00057 #include "input.h"
00058 #include "taulines.h"
00059 #include "version.h"
00060 #include "thermal.h"
00061 #include "optimize.h"
00062 #include "grid.h"
00063 #include "timesc.h"
00064 #include "cloudy.h"
00065 #include "warnings.h"
00066 #include "lines_service.h"
00067 #include "cddrive.h"
00068
00069
00070
00071
00072
00073
00074
00075 int cdDrive()
00076 {
00077 bool lgBAD;
00078
00079 DEBUG_ENTRY( "cdDrive()" );
00080
00081
00082
00083
00084
00085
00086 if( !lgcdInitCalled )
00087 {
00088 printf(" cdInit was not called first - this must be the first call.\n");
00089 cdEXIT(EXIT_FAILURE);
00090 }
00091
00092 if( trace.lgTrace )
00093 {
00094 fprintf( ioQQQ,
00095 "cdDrive: lgOptimr=%1i lgVaryOn=%1i lgNoVary=%1i input.nSave:%li\n",
00096 optimize.lgOptimr , optimize.lgVaryOn , optimize.lgNoVary, input.nSave );
00097 }
00098
00099
00100
00101
00102 if( optimize.lgOptimr && optimize.lgVaryOn && !optimize.lgNoVary )
00103 optimize.lgVaryOn = true;
00104 else
00105 optimize.lgVaryOn = false;
00106
00107
00108
00109
00110 InitCoreload();
00111
00112 if( optimize.lgVaryOn )
00113 {
00114
00115 if( trace.lgTrace )
00116 fprintf( ioQQQ, "cdDrive: calling grid_do\n");
00117
00118 lgBAD = grid_do();
00119 }
00120 else
00121 {
00122 if( trace.lgTrace )
00123 fprintf( ioQQQ, "cdDrive: calling cloudy\n");
00124
00125
00126 lgBAD = cloudy();
00127 }
00128
00129
00130 lgcdInitCalled = false;
00131
00132 if( lgAbort || lgBAD )
00133 {
00134 if( trace.lgTrace )
00135 fprintf( ioQQQ, "cdDrive: returning failure during call. \n");
00136
00137 return 1;
00138 }
00139 else
00140 {
00141
00142 return 0;
00143 }
00144 }
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 void cdPrtWL( FILE *io , realnum wl )
00156 {
00157
00158 DEBUG_ENTRY( "cdPrtWL()" );
00159
00160 prt_wl( io , wl );
00161 return;
00162 }
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173 void cdReasonGeo(FILE * ioOUT)
00174 {
00175
00176 DEBUG_ENTRY( "cdReasonGeo()" );
00177
00178
00179 fprintf( ioOUT, "%s", warnings.chRgcln[0] );
00180 fprintf( ioOUT , "\n" );
00181
00182 fprintf( ioOUT, "%s", warnings.chRgcln[1] );
00183 fprintf( ioOUT , "\n" );
00184 return;
00185 }
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196 void cdWarnings(FILE *ioPNT )
00197 {
00198 long int i;
00199
00200 DEBUG_ENTRY( "cdWarnings()" );
00201
00202 if( warnings.nwarn > 0 )
00203 {
00204 for( i=0; i < warnings.nwarn; i++ )
00205 {
00206
00207 fprintf( ioPNT, "%s", warnings.chWarnln[i] );
00208 fprintf( ioPNT, "\n" );
00209 }
00210 }
00211
00212 return;
00213 }
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 void cdCautions(FILE * ioOUT)
00225 {
00226 long int i;
00227
00228 DEBUG_ENTRY( "cdCautions()" );
00229
00230 if( warnings.ncaun > 0 )
00231 {
00232 for( i=0; i < warnings.ncaun; i++ )
00233 {
00234 fprintf( ioOUT, "%s", warnings.chCaunln[i] );
00235 fprintf( ioOUT, "\n" );
00236 }
00237 }
00238 return;
00239 }
00240
00241
00242
00243
00244
00245
00246
00247 void cdTimescales(
00248
00249 double *TTherm ,
00250
00251 double *THRecom ,
00252
00253 double *TH2 )
00254 {
00255
00256 DEBUG_ENTRY( "cdTimescales()" );
00257
00258
00259
00260
00261 *TTherm = timesc.time_therm_long;
00262
00263
00264 *THRecom = timesc.time_Hrecom_long;
00265
00266
00267 *TH2 = MAX2( timesc.time_H2_Dest_longest , timesc.time_H2_Form_longest );
00268 return;
00269 }
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280 void cdSurprises(FILE * ioOUT)
00281 {
00282 long int i;
00283
00284 DEBUG_ENTRY( "cdSurprises()" );
00285
00286 if( warnings.nbang > 0 )
00287 {
00288 for( i=0; i < warnings.nbang; i++ )
00289 {
00290 fprintf( ioOUT, "%s", warnings.chBangln[i] );
00291 fprintf( ioOUT, "\n" );
00292 }
00293 }
00294
00295 return;
00296 }
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307 void cdNotes(FILE * ioOUT)
00308 {
00309 long int i;
00310
00311 DEBUG_ENTRY( "cdNotes()" );
00312
00313 if( warnings.nnote > 0 )
00314 {
00315 for( i=0; i < warnings.nnote; i++ )
00316 {
00317 fprintf( ioOUT, "%s", warnings.chNoteln[i] );
00318 fprintf( ioOUT, "\n" );
00319 }
00320 }
00321 return;
00322 }
00323
00324
00325
00326
00327
00328
00329
00330
00331 double cdCooling_last()
00332 {
00333 return thermal.ctot;
00334 }
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344 void cdVersion(char chString[])
00345 {
00346 strcpy( chString , t_version::Inst().chVersion );
00347 return;
00348 }
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359 void cdDate(char chString[])
00360 {
00361 strcpy( chString , t_version::Inst().chDate );
00362 return;
00363 }
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374 double cdHeating_last()
00375 {
00376 return thermal.htot;
00377 }
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388 double cdEDEN_last()
00389 {
00390 return dense.eden;
00391 }
00392
00393
00394
00395
00396
00397
00398
00399
00400 #include "noexec.h"
00401
00402 void cdNoExec()
00403 {
00404
00405 DEBUG_ENTRY( "cdNoExec()" );
00406
00407
00408
00409 noexec.lgNoExec = true;
00410
00411 return;
00412 }
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423 static bool lgCalled=false;
00424
00425
00426 #if defined(_MSC_VER)
00427
00428
00429 struct timeval {
00430 long tv_sec;
00431 long tv_usec;
00432 };
00433 #else
00434 #include <sys/time.h>
00435 #include <sys/resource.h>
00436 #endif
00437
00438
00439 static struct timeval before;
00440
00441
00442 STATIC void cdClock(struct timeval *clock_dat)
00443 {
00444 DEBUG_ENTRY( "cdClock()" );
00445
00446
00447
00448 #if defined(_MSC_VER) || defined(__HP_aCC)
00449 clock_t clock_val;
00450
00451
00452 clock_val = clock();
00453 clock_dat->tv_sec = clock_val/CLOCKS_PER_SEC;
00454 clock_dat->tv_usec = 1000000*(clock_val-(clock_dat->tv_sec*CLOCKS_PER_SEC))/CLOCKS_PER_SEC;
00455
00456
00457 #else
00458 struct rusage rusage;
00459 if(getrusage(RUSAGE_SELF,&rusage) != 0)
00460 {
00461 fprintf( ioQQQ, "DISASTER cdClock called getrusage with invalid arguments.\n" );
00462 fprintf( ioQQQ, "Sorry.\n" );
00463 cdEXIT(EXIT_FAILURE);
00464 }
00465 clock_dat->tv_sec = rusage.ru_utime.tv_sec;
00466 clock_dat->tv_usec = rusage.ru_utime.tv_usec;
00467 #endif
00468
00469 }
00470
00471
00472 void cdSetExecTime()
00473 {
00474 cdClock(&before);
00475 lgCalled = true;
00476 }
00477
00478
00479 double cdExecTime()
00480 {
00481 DEBUG_ENTRY( "cdExecTime()" );
00482
00483 struct timeval clock_dat;
00484
00485 if( lgCalled)
00486 {
00487 cdClock(&clock_dat);
00488
00489
00490
00491 return (double)(clock_dat.tv_sec-before.tv_sec)+1e-6*(double)(clock_dat.tv_usec-before.tv_usec);
00492 }
00493 else
00494 {
00495
00496
00497 fprintf( ioQQQ, "DISASTER cdExecTime was called before SetExecTime, impossible.\n" );
00498 fprintf( ioQQQ, "Sorry.\n" );
00499 cdEXIT(EXIT_FAILURE);
00500 }
00501 }
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511 void cdPrintCommands( FILE * ioOUT )
00512 {
00513 long int i;
00514 fprintf( ioOUT, " Input commands follow:\n" );
00515 fprintf( ioOUT, "c ======================\n" );
00516
00517 for( i=0; i <= input.nSave; i++ )
00518 {
00519 fprintf( ioOUT, "%s\n", input.chCardSav[i] );
00520 }
00521 fprintf( ioOUT, "c ======================\n" );
00522 }
00523
00524
00525
00526
00527
00528
00529
00530
00531 long int cdEmis(
00532
00533
00534
00535 char *chLabel,
00536
00537 realnum wavelength,
00538
00539 double *emiss )
00540 {
00541 DEBUG_ENTRY( "cdEms()" );
00542 long int i = cdEmis( chLabel , wavelength , emiss, false );
00543 return i;
00544 }
00545
00546
00547
00548 long int cdEmis(
00549
00550
00551
00552 char *chLabel,
00553
00554 realnum wavelength,
00555
00556 double *emiss ,
00557
00558 bool lgEmergent )
00559 {
00560
00561 char chCARD[INPUT_LINE_LENGTH];
00562 char chCaps[5];
00563 long int j;
00564 realnum errorwave;
00565
00566 DEBUG_ENTRY( "cdEms()" );
00567
00568
00569
00570
00571
00572 long int ipEmisType = 0;
00573 if( lgEmergent )
00574 ipEmisType = 1;
00575
00576 strcpy( chCARD, chLabel );
00577
00578
00579 caps(chCARD);
00580
00581
00582 errorwave = WavlenErrorGet( wavelength );
00583
00584 for( j=0; j < LineSave.nsum; j++ )
00585 {
00586
00587 cap4(chCaps, LineSv[j].chALab);
00588
00589
00590
00591
00592 if( fabs(LineSv[j].wavelength-wavelength) < errorwave && strcmp(chCaps,chCARD) == 0 )
00593 {
00594
00595 *emiss = LineSv[j].emslin[ipEmisType];
00596
00597 return j;
00598 }
00599 }
00600
00601
00602 return -LineSave.nsum;
00603 }
00604
00605
00606
00607
00608
00609
00610
00611 void cdEmis_ip(
00612
00613 long int ipLine,
00614
00615 double *emiss ,
00616
00617 bool lgEmergent )
00618 {
00619 DEBUG_ENTRY( "cdEmis_ip()" );
00620
00621
00622
00623 ASSERT( ipLine >= 0 && ipLine < LineSave.nsum );
00624 *emiss = LineSv[ipLine].emslin[lgEmergent];
00625 return;
00626 }
00627
00628
00629
00630
00631
00632
00633
00634 int cdColm(
00635
00636
00637
00638 const char *chLabel,
00639
00640
00641
00642 long int ion,
00643
00644
00645 double *theocl )
00646 {
00647 long int nelem;
00648
00649 char chLABEL_CAPS[20];
00650
00651 DEBUG_ENTRY( "cdColm()" );
00652
00653
00654 if( chLabel[4] != '\0' || strlen(chLabel) != 4 )
00655 {
00656 fprintf( ioQQQ, " cdColm called with insane chLabel (between quotes) \"%s\", must be exactly 4 characters long.\n",
00657 chLabel );
00658 return 1;
00659 }
00660
00661 strcpy( chLABEL_CAPS, chLabel );
00662
00663 caps(chLABEL_CAPS);
00664
00665
00666
00667
00668 if( ion < 0 )
00669 {
00670 fprintf( ioQQQ, " cdColm called with insane ion, =%li\n",
00671 ion );
00672 return 1;
00673 }
00674
00675 else if( ion == 0 )
00676 {
00677 struct molecule *sp;
00678
00679
00680 if( strcmp( chLABEL_CAPS , "H2 " )==0 )
00681 {
00682 *theocl = colden.colden[ipCOL_H2g] + colden.colden[ipCOL_H2s];
00683 }
00684
00685
00686 else if( strcmp( chLABEL_CAPS , "H- " )==0 )
00687 {
00688 *theocl = colden.colden[ipCOL_HMIN];
00689 }
00690
00691
00692 else if( strcmp( chLABEL_CAPS , "H2+ " )==0 )
00693 {
00694 *theocl = colden.colden[ipCOL_H2p];
00695 }
00696
00697
00698 else if( strcmp( chLABEL_CAPS , "H3+ " )==0 )
00699 {
00700 *theocl = colden.colden[ipCOL_H3p];
00701 }
00702
00703
00704 else if( strcmp( chLABEL_CAPS , "H2G " )==0 )
00705 {
00706 *theocl = colden.colden[ipCOL_H2g];
00707 }
00708
00709
00710 else if( strcmp( chLABEL_CAPS , "H2* " )==0 )
00711 {
00712 *theocl = colden.colden[ipCOL_H2s];
00713 }
00714
00715
00716 else if( strcmp( chLABEL_CAPS , "HEH+" )==0 )
00717 {
00718 *theocl = colden.colden[ipCOL_HeHp];
00719 }
00720
00721
00722 else if( (sp = findspecies(chLabel)) != &null_mole )
00723 {
00724 *theocl = sp->hevcol;
00725 }
00726
00727
00728
00729
00730
00731 else if( strcmp( chLABEL_CAPS , "CII*" )==0 )
00732 {
00733 *theocl = colden.C2Colden[1];
00734 }
00735 else if( strcmp( chLABEL_CAPS , "C11*" )==0 )
00736 {
00737 *theocl = colden.C1Colden[0];
00738 }
00739 else if( strcmp( chLABEL_CAPS , "C12*" )==0 )
00740 {
00741 *theocl = colden.C1Colden[1];
00742 }
00743 else if( strcmp( chLABEL_CAPS , "C13*" )==0 )
00744 {
00745 *theocl = colden.C1Colden[2];
00746 }
00747 else if( strcmp( chLABEL_CAPS , "O11*" )==0 )
00748 {
00749 *theocl = colden.O1Colden[0];
00750 }
00751 else if( strcmp( chLABEL_CAPS , "O12*" )==0 )
00752 {
00753 *theocl = colden.O1Colden[1];
00754 }
00755 else if( strcmp( chLABEL_CAPS , "O13*" )==0 )
00756 {
00757 *theocl = colden.O1Colden[2];
00758 }
00759
00760 else if( strcmp( chLABEL_CAPS , "C30*" )==0 )
00761 {
00762 *theocl = colden.C3Colden[1];
00763 }
00764 else if( strcmp( chLABEL_CAPS , "C31*" )==0 )
00765 {
00766 *theocl = colden.C3Colden[2];
00767 }
00768 else if( strcmp( chLABEL_CAPS , "C32*" )==0 )
00769 {
00770 *theocl = colden.C3Colden[3];
00771 }
00772 else if( strcmp( chLABEL_CAPS , "SI2*" )==0 )
00773 {
00774 *theocl = colden.Si2Colden[1];
00775 }
00776 else if( strcmp( chLABEL_CAPS , "HE1*" )==0 )
00777 {
00778 *theocl = colden.He123S;
00779 }
00780
00781 else if( strncmp(chLABEL_CAPS , "H2" , 2 ) == 0 )
00782 {
00783 long int iVib = chLABEL_CAPS[2] - '0';
00784 long int iRot = chLABEL_CAPS[3] - '0';
00785 if( iVib<0 || iRot < 0 )
00786 {
00787 fprintf( ioQQQ, " cdColm called with insane v,J for H2=\"%4.4s\" caps=\"%4.4s\"\n",
00788 chLabel , chLABEL_CAPS );
00789 return 1;
00790 }
00791 *theocl = cdH2_colden( iVib , iRot );
00792 }
00793
00794
00795 else
00796 {
00797 fprintf( ioQQQ, " cdColm called with unknown element chLabel, org=\"%4.4s \" 0 caps=\"%4.4s\" 0\n",
00798 chLabel , chLABEL_CAPS );
00799 return 1;
00800 }
00801 }
00802 else
00803 {
00804
00805
00806 nelem = 0;
00807 while( nelem < LIMELM &&
00808 strncmp(chLABEL_CAPS,elementnames.chElementNameShort[nelem],4) != 0 )
00809 {
00810 ++nelem;
00811 }
00812
00813
00814
00815 if( nelem < LIMELM )
00816 {
00817
00818
00819
00820 if( ion > MAX2(3,nelem + 2) )
00821 {
00822 fprintf( ioQQQ,
00823 " cdColm asked to return ionization stage %ld for element %s but this is not physical.\n",
00824 ion, chLabel );
00825 return 1;
00826 }
00827
00828
00829 *theocl = mean.xIonMean[0][nelem][ion-1][0];
00830
00831
00832
00833
00834 if( nelem==ipHYDROGEN && ion==3 )
00835 *theocl /= 2.;
00836 }
00837 else
00838 {
00839 fprintf( ioQQQ,
00840 " cdColm did not understand this combination of ion %4ld and element %4.4s.\n",
00841 ion, chLabel );
00842 return 1;
00843 }
00844 }
00845 return 0;
00846 }
00847
00848
00849
00850
00851
00852
00853
00854
00855 void cdErrors(FILE *ioOUT)
00856 {
00857 long int nc,
00858 nn,
00859 npe,
00860 ns,
00861 nte,
00862 nw ,
00863 nIone,
00864 nEdene;
00865 bool lgAbort_loc;
00866
00867 DEBUG_ENTRY( "cdErrors()" );
00868
00869
00870 cdNwcns(&lgAbort_loc,&nw,&nc,&nn,&ns,&nte,&npe, &nIone, &nEdene );
00871
00872
00873 if( nw || nc || nte || npe || nIone || nEdene || lgAbort_loc )
00874 {
00875
00876 fprintf( ioOUT, "%75.75s\n", input.chTitle );
00877
00878 if( lgAbort_loc )
00879 fprintf(ioOUT," Calculation ended with abort!\n");
00880
00881
00882 if( nw != 0 )
00883 {
00884 cdWarnings(ioOUT);
00885 }
00886
00887
00888 if( nc != 0 )
00889 {
00890 cdCautions(ioOUT);
00891 }
00892
00893 if( nte != 0 )
00894 {
00895 fprintf( ioOUT , "Te failures=%4ld\n", nte );
00896 }
00897
00898 if( npe != 0 )
00899 {
00900 fprintf( ioOUT , "Pressure failures=%4ld\n", npe );
00901 }
00902
00903 if( nIone != 0 )
00904 {
00905 fprintf( ioOUT , "Ionization failures=%4ld\n", nte );
00906 }
00907
00908 if( nEdene != 0 )
00909 {
00910 fprintf( ioOUT , "Electron density failures=%4ld\n", npe );
00911 }
00912 }
00913 return;
00914 }
00915
00916
00917
00918
00919
00920
00921 void cdDepth_depth( double cdDepth[] )
00922 {
00923 long int nz;
00924
00925 DEBUG_ENTRY( "cdDepth_depth()" );
00926
00927 for( nz = 0; nz<nzone; ++nz )
00928 {
00929 cdDepth[nz] = struc.depth[nz];
00930 }
00931 return;
00932 }
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946 void cdPressure_depth(
00947
00948 double TotalPressure[],
00949
00950 double GasPressure[],
00951
00952 double RadiationPressure[])
00953 {
00954 long int nz;
00955
00956 DEBUG_ENTRY( "cdPressure_depth()" );
00957
00958 for( nz = 0; nz<nzone; ++nz )
00959 {
00960 TotalPressure[nz] = struc.pressure[nz];
00961 GasPressure[nz] = struc.GasPressure[nz];
00962 RadiationPressure[nz] = struc.pres_radiation_lines_curr[nz];
00963 }
00964 return;
00965 }
00966
00967
00968
00969
00970
00971
00972
00973 void cdPressure_last(
00974 double *PresTotal,
00975 double *PresGas,
00976 double *PresRad)
00977 {
00978
00979 DEBUG_ENTRY( "cdPressure_last()" );
00980
00981 *PresGas = pressure.PresGasCurr;
00982 *PresRad = pressure.pres_radiation_lines_curr;
00983 *PresTotal = pressure.PresTotlCurr;
00984 return;
00985 }
00986
00987
00988
00989
00990
00991
00992
00993
00994 long int cdnZone()
00995 {
00996 return nzone;
00997 }
00998
00999
01000
01001
01002
01003
01004
01005
01006 double cdTemp_last()
01007 {
01008 return phycon.te;
01009 }
01010
01011
01012
01013
01014
01015
01016
01017 int cdIonFrac(
01018
01019 const char *chLabel,
01020
01021
01022 long int IonStage,
01023
01024 double *fracin,
01025
01026 const char *chWeight ,
01027
01028 bool lgDensity )
01029
01030 {
01031 long int ip,
01032 ion,
01033 nelem;
01034 realnum aaa[LIMELM + 1];
01035
01036 char chCARD[INPUT_LINE_LENGTH];
01037
01038 DEBUG_ENTRY( "cdIonFrac()" );
01039
01040 strcpy( chCARD, chWeight );
01041
01042 caps(chCARD);
01043
01044 int dim;
01045 if( strcmp(chCARD,"RADIUS") == 0 )
01046 dim = 0;
01047 else if( strcmp(chCARD,"AREA") == 0 )
01048 dim = 1;
01049 else if( strcmp(chCARD,"VOLUME") == 0 )
01050 dim = 2;
01051 else
01052 {
01053 fprintf( ioQQQ, " cdIonFrac: chWeight=%6.6s makes no sense to me, valid options are RADIUS, AREA, and VOLUME\n",
01054 chWeight );
01055 *fracin = 0.;
01056 return 1;
01057 }
01058
01059
01060 strcpy( chCARD, chLabel );
01061
01062 caps(chCARD);
01063
01064 if( IonStage==0 )
01065 {
01066
01067 if( strcmp(chCARD,"H2 " ) == 0 )
01068 {
01069
01070 nelem = 0;
01071 IonStage = 3;
01072 }
01073 else
01074 {
01075 fprintf( ioQQQ, " cdIonFrac: ion stage of zero and element %s makes no sense to me\n",
01076 chCARD );
01077 *fracin = 0.;
01078 return 1;
01079 }
01080 }
01081
01082 else
01083 {
01084
01085 nelem = 0;
01086 while( nelem < LIMELM &&
01087 strcmp(chCARD,elementnames.chElementNameShort[nelem]) != 0 )
01088 {
01089 ++nelem;
01090 }
01091 }
01092
01093
01094 if( nelem >= LIMELM )
01095 {
01096 fprintf( ioQQQ, " cdIonFrac called with unknown element chLabel, =%4.4s\n",
01097 chLabel );
01098 return 1;
01099 }
01100
01101
01102
01103
01104
01105
01106 ion = IonStage - 1;
01107
01108 if( (ion > nelem+1 || ion < 0 ) && !(nelem==ipHYDROGEN&&ion==2))
01109 {
01110 fprintf( ioQQQ, " cdIonFrac asked to return ionization stage %ld for element %4.4s but this is not physical.\n",
01111 IonStage, chLabel );
01112 *fracin = -1.;
01113 return 1;
01114 }
01115
01116
01117
01118
01119
01120
01121 mean.MeanIon('i',nelem,dim,&ip,aaa,lgDensity);
01122 *fracin = pow((realnum)10.f,aaa[ion]);
01123
01124
01125 return 0;
01126 }
01127
01128
01129
01130
01131
01132
01133
01134
01135
01136
01137
01138
01139 long debugLine( realnum wavelength )
01140 {
01141 long j, kount;
01142 realnum errorwave;
01143
01144 kount = 0;
01145
01146
01147 errorwave = WavlenErrorGet( wavelength );
01148
01149 for( j=0; j < LineSave.nsum; j++ )
01150 {
01151
01152
01153 if( fabs(LineSv[j].wavelength-wavelength) < errorwave )
01154 {
01155 printf("%s\n", LineSv[j].chALab);
01156 ++kount;
01157 }
01158 }
01159 printf(" hits = %li\n", kount );
01160 return kount;
01161 }
01162
01163
01164
01165
01166
01167
01168
01169
01170
01171
01172 long int cdLine(
01173 const char *chLabel,
01174
01175 realnum wavelength,
01176
01177 double *relint,
01178
01179 double *absint )
01180 {
01181 DEBUG_ENTRY( "cdLine()" );
01182 long int i = cdLine( chLabel , wavelength , relint , absint, false );
01183 return i;
01184 }
01185 long int cdLine(
01186 const char *chLabel,
01187
01188 realnum wavelength,
01189
01190 double *relint,
01191
01192 double *absint ,
01193 bool lgEmergent )
01194 {
01195 char chCaps[5],
01196 chFind[5];
01197 long int ipobs,
01198 j, index_of_closest=LONG_MIN,
01199 index_of_closest_w_correct_label=-1;
01200 realnum errorwave, smallest_error=BIGFLOAT,
01201 smallest_error_w_correct_label=BIGFLOAT;
01202
01203 DEBUG_ENTRY( "cdLine()" );
01204
01205
01206 if( LineSave.nsum == 0 )
01207 {
01208 *relint = 0.;
01209 *absint = 0.;
01210 return 0;
01211 }
01212 ASSERT( LineSave.ipNormWavL >= 0);
01213 ASSERT( LineSave.nsum > 0);
01214
01215
01216 if( chLabel[4] != '\0' || strlen(chLabel) != 4 )
01217 {
01218 fprintf( ioQQQ, " cdLine called with insane chLabel (between quotes) \"%s\", must be exactly 4 characters long.\n",
01219 chLabel );
01220 return 1;
01221 }
01222
01223
01224 cap4(chFind, chLabel);
01225
01226
01227
01228 for( j=0; j<=3; j++ )
01229 {
01230 if( chFind[j] == '\t' )
01231 {
01232 chFind[j] = ' ';
01233 }
01234 }
01235
01236
01237 errorwave = WavlenErrorGet( wavelength );
01238
01239 {
01240
01241 enum{DEBUG_LOC=false};
01242
01243 if( DEBUG_LOC && fabs(wavelength-1000.)<0.01 )
01244 {
01245 fprintf(ioQQQ,"cdDrive wl %.4e error %.3e\n",
01246 wavelength, errorwave );
01247 }
01248 }
01249
01250
01251 for( j=1; j < LineSave.nsum; j++ )
01252 {
01253
01254
01255 realnum current_error;
01256 current_error = (realnum)fabs(LineSv[j].wavelength-wavelength);
01257
01258
01259 cap4(chCaps, LineSv[j].chALab);
01260
01261 if( current_error < smallest_error )
01262 {
01263 index_of_closest = j;
01264 smallest_error = current_error;
01265 }
01266
01267 if( current_error < smallest_error_w_correct_label &&
01268 (strcmp(chCaps,chFind) == 0) )
01269 {
01270 index_of_closest_w_correct_label = j;
01271 smallest_error_w_correct_label = current_error;
01272 }
01273
01274
01275
01276 if( current_error <= errorwave ||
01277 fp_equal( wavelength + errorwave, LineSv[j].wavelength ) ||
01278 fp_equal( wavelength - errorwave, LineSv[j].wavelength ) )
01279 {
01280
01281 if( strcmp(chCaps,chFind) == 0 )
01282 {
01283
01284 ipobs = j;
01285
01286
01287 if( LineSv[LineSave.ipNormWavL].SumLine[lgEmergent] > 0. )
01288 {
01289 *relint = LineSv[ipobs].SumLine[lgEmergent]/
01290 LineSv[LineSave.ipNormWavL].SumLine[lgEmergent]*
01291 LineSave.ScaleNormLine;
01292 }
01293 else
01294 {
01295 *relint = 0.;
01296 }
01297
01298
01299 if( LineSv[ipobs].SumLine[lgEmergent] > 0. )
01300 {
01301 *absint = log10(LineSv[ipobs].SumLine[lgEmergent]) +
01302 radius.Conv2PrtInten;
01303 }
01304 else
01305 {
01306
01307 *absint = -37.;
01308 }
01309
01310 return ipobs;
01311 }
01312 }
01313 }
01314
01315
01316
01317 fprintf( ioQQQ," PROBLEM cdLine did not find line with label (between quotes) \"%4s\" and wavelength ", chFind );
01318 prt_wl(ioQQQ,wavelength);
01319 if( index_of_closest >= 0 )
01320 {
01321 fprintf( ioQQQ,".\n The closest line (any label) was \"%4s\"\t",
01322 LineSv[index_of_closest].chALab );
01323 prt_wl(ioQQQ,LineSv[index_of_closest].wavelength );
01324 if( index_of_closest_w_correct_label >= 0 )
01325 {
01326 fprintf( ioQQQ,"\n The closest with correct label was \"%4s\"\t", chFind );
01327 prt_wl(ioQQQ,LineSv[index_of_closest_w_correct_label].wavelength );
01328 fprintf( ioQQQ,"\n" );
01329 }
01330 else
01331 fprintf( ioQQQ,"\n No line found with label \"%s\".\n", chFind );
01332 fprintf( ioQQQ,"\n" );
01333 }
01334 else
01335 {
01336 fprintf( ioQQQ,".\n PROBLEM No close line was found\n" );
01337 TotalInsanity();
01338 }
01339
01340 *absint = 0.;
01341 *relint = 0.;
01342
01343
01344
01345 return -LineSave.nsum;
01346 }
01347
01348
01349 void cdLine_ip(long int ipLine,
01350
01351 double *relint,
01352
01353 double *absint )
01354 {
01355
01356 DEBUG_ENTRY( "cdLine_ip()" );
01357 cdLine_ip( ipLine , relint , absint , false );
01358 }
01359 void cdLine_ip(long int ipLine,
01360
01361 double *relint,
01362
01363 double *absint ,
01364 bool lgEmergent )
01365 {
01366
01367 DEBUG_ENTRY( "cdLine_ip()" );
01368
01369
01370 if( LineSave.nsum == 0 )
01371 {
01372 *relint = 0.;
01373 *absint = 0.;
01374 return;
01375 }
01376 ASSERT( LineSave.ipNormWavL >= 0);
01377 ASSERT( LineSave.nsum > 0);
01378
01379
01380 if( LineSv[LineSave.ipNormWavL].SumLine[lgEmergent] > 0. )
01381 *relint = LineSv[ipLine].SumLine[lgEmergent]/
01382 LineSv[LineSave.ipNormWavL].SumLine[lgEmergent]*
01383 LineSave.ScaleNormLine;
01384 else
01385 *relint = 0.;
01386
01387
01388 if( LineSv[ipLine].SumLine[lgEmergent] > 0. )
01389 *absint = log10(LineSv[ipLine].SumLine[lgEmergent]) +
01390 radius.Conv2PrtInten;
01391 else
01392
01393 *absint = -37.;
01394
01395 return;
01396 }
01397
01398
01399
01400
01401
01402
01403
01404 void cdNwcns(
01405
01406 bool *lgAbort_ret ,
01407
01408 long int *NumberWarnings,
01409 long int *NumberCautions,
01410 long int *NumberNotes,
01411 long int *NumberSurprises,
01412
01413 long int *NumberTempFailures,
01414
01415 long int *NumberPresFailures,
01416
01417 long int *NumberIonFailures,
01418
01419 long int *NumberNeFailures )
01420 {
01421
01422 DEBUG_ENTRY( "cdNwcns()" );
01423
01424
01425 *lgAbort_ret = lgAbort;
01426
01427 *NumberWarnings = warnings.nwarn;
01428 *NumberCautions = warnings.ncaun;
01429 *NumberNotes = warnings.nnote;
01430 *NumberSurprises = warnings.nbang;
01431
01432
01433 *NumberTempFailures = conv.nTeFail;
01434 *NumberPresFailures = conv.nPreFail;
01435 *NumberIonFailures = conv.nIonFail;
01436 *NumberNeFailures = conv.nNeFail;
01437 return;
01438 }
01439
01440 void cdOutput( const char* filename, const char *mode )
01441 {
01442 DEBUG_ENTRY( "cdOutput()" );
01443
01444 if( ioQQQ != stdout && ioQQQ != NULL )
01445 fclose(ioQQQ);
01446 FILE *fp = stdout;
01447 if( *filename != '\0' )
01448 fp = open_data( filename, mode, AS_LOCAL_ONLY );
01449 ioQQQ = fp;
01450 }
01451
01452 void cdInput( const char* filename, const char *mode )
01453 {
01454 DEBUG_ENTRY( "cdInput()" );
01455
01456 if( ioStdin != stdin && ioStdin != NULL )
01457 fclose(ioStdin);
01458 FILE *fp = stdin;
01459 if( *filename != '\0' )
01460 fp = open_data( filename, mode, AS_LOCAL_ONLY );
01461 ioStdin = fp;
01462 }
01463
01464
01465
01466
01467
01468
01469
01470 void cdTalk(bool lgTOn)
01471 {
01472
01473 DEBUG_ENTRY( "cdTalk()" );
01474
01475
01476 called.lgTalk = lgTOn && cpu.lgMPI_talk();
01477
01478 called.lgTalkForcedOff = !lgTOn;
01479 return;
01480 }
01481
01482
01483 double cdB21cm()
01484 {
01485 double ret;
01486
01487 DEBUG_ENTRY( "cdB21cm()" );
01488
01489
01490 if( mean.TempB_HarMean[0][1] > SMALLFLOAT )
01491 {
01492 ret = mean.TempB_HarMean[0][0]/mean.TempB_HarMean[0][1];
01493 }
01494 else
01495 {
01496 ret = 0.;
01497 }
01498 return ret;
01499 }
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516
01517
01518
01522
01523
01524 int cdTemp(
01525
01526 const char *chLabel,
01527
01528
01529 long int IonStage,
01530
01531 double *TeMean,
01532
01533 const char *chWeight )
01534 {
01535 long int ip,
01536 ion,
01537 nelem;
01538 realnum aaa[LIMELM + 1];
01539
01540 char chWGHT[INPUT_LINE_LENGTH] , chELEM[INPUT_LINE_LENGTH];
01541
01542 DEBUG_ENTRY( "cdTemp()" );
01543
01544
01545 strcpy( chWGHT, chWeight );
01546 caps(chWGHT);
01547 strcpy( chELEM, chLabel );
01548 caps(chELEM);
01549
01550
01551 int dim;
01552 if( strcmp(chWGHT,"RADIUS") == 0 )
01553 dim = 0;
01554 else if( strcmp(chWGHT,"AREA") == 0 )
01555 dim = 1;
01556 else if( strcmp(chWGHT,"VOLUME") == 0 )
01557 dim = 2;
01558 else
01559 {
01560 fprintf( ioQQQ, " cdTemp: chWeight=%6.6s makes no sense to me, the options are RADIUS, AREA, and VOLUME.\n",
01561 chWeight );
01562 *TeMean = 0.;
01563 return 1;
01564 }
01565
01566 if( IonStage == 0 )
01567 {
01568
01569 if( strcmp(chELEM,"21CM") == 0 )
01570 {
01571 if( mean.TempHarMean[dim][1] > SMALLFLOAT )
01572 *TeMean = mean.TempHarMean[dim][0]/mean.TempHarMean[dim][1];
01573 else
01574 *TeMean = 0.;
01575 }
01576
01577
01578 else if( strcmp(chELEM,"SPIN") == 0 )
01579 {
01580 *TeMean = mean.TempH_21cmSpinMean[dim][0] / SDIV(mean.TempH_21cmSpinMean[dim][1]);
01581 }
01582
01583 else if( strcmp(chELEM,"OPTI") == 0 )
01584 {
01585 *TeMean =
01586 3.84e-7 * Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->TauCon /
01587 SDIV( HFLines[0].Emis->TauCon );
01588 }
01589
01590 else if( strcmp(chELEM,"H2 ") == 0 )
01591 {
01592 if( mean.TempIonMean[dim][ipHYDROGEN][2][1] > SMALLFLOAT )
01593 *TeMean = mean.TempIonMean[dim][ipHYDROGEN][2][0] /
01594 mean.TempIonMean[dim][ipHYDROGEN][2][1];
01595
01596 else
01597 *TeMean = 0.;
01598 }
01599
01600 else if( strcmp(chELEM,"TENE") == 0 )
01601 {
01602 if( mean.TempEdenMean[dim][1] > SMALLFLOAT )
01603 *TeMean = mean.TempEdenMean[dim][0]/mean.TempEdenMean[dim][1];
01604 else
01605 *TeMean = 0.;
01606 }
01607
01608 else if( strcmp(chELEM," ") == 0 )
01609 {
01610 if( mean.TempMean[dim][1] > SMALLFLOAT )
01611 *TeMean = mean.TempMean[dim][0]/mean.TempMean[dim][1];
01612 else
01613 *TeMean = 0.;
01614 }
01615 else
01616 {
01617 fprintf( ioQQQ, " cdTemp called with ion=0 and unknown quantity, =%4.4s\n",
01618 chLabel );
01619 return 1;
01620 }
01621
01622
01623 return 0;
01624 }
01625
01626
01627 nelem = 0;
01628 while( nelem < LIMELM &&
01629 strcmp(chELEM,elementnames.chElementNameShort[nelem]) != 0 )
01630 {
01631 ++nelem;
01632 }
01633
01634
01635 if( nelem >= LIMELM )
01636 {
01637 fprintf( ioQQQ, " cdTemp called with unknown element chLabel, =%4.4s\n",
01638 chLabel );
01639 return 1;
01640 }
01641
01642
01643
01644
01645
01646
01647 ion = IonStage - 1;
01648
01649 if( ion > nelem+1 || ion < 0 )
01650 {
01651 fprintf( ioQQQ, " cdTemp asked to return ionization stage %ld for element %4.4s but this is not physical.\n",
01652 IonStage, chLabel );
01653 return 1;
01654 }
01655
01656
01657
01658
01659 mean.MeanIon('t', nelem,dim,&ip,aaa,false);
01660 *TeMean = pow((realnum)10.f,aaa[ion]);
01661 return 0;
01662 }
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672
01673
01674 int cdRead(
01675
01676 const char *chInputLine )
01677 {
01678 char *chEOL ,
01679 chCARD[INPUT_LINE_LENGTH],
01680 chLocal[INPUT_LINE_LENGTH];
01681 bool lgComment;
01682
01683 DEBUG_ENTRY( "cdRead()" );
01684
01685
01686
01687 if( !lgcdInitCalled )
01688 {
01689 printf(" cdInit was not called first - this must be the first call.\n");
01690 cdEXIT(EXIT_FAILURE);
01691 }
01692
01693
01694
01695
01696 if( ( lgInputComment( chInputLine ) ||
01697
01698 chInputLine[0] == '\n' || chInputLine[0] == ' ' ) &&
01699
01700 ! ( chInputLine[0] == 'C' || chInputLine[0] == 'c' ) )
01701 {
01702
01703 return NKRD - input.nSave;
01704 }
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714 ++input.nSave;
01715
01716 if( input.nSave >= NKRD )
01717 {
01718
01719 fprintf( ioQQQ,
01720 " Too many line images entered to cdRead. The limit is %d\n",
01721 NKRD );
01722 fprintf( ioQQQ,
01723 " The limit to the number of allowed input lines is %d. This limit was exceeded. Sorry.\n",
01724 NKRD );
01725 fprintf( ioQQQ,
01726 " This limit is set by the variable NKRD which appears in input.h \n" );
01727 fprintf( ioQQQ,
01728 " Increase it everywhere it appears.\n" );
01729 cdEXIT(EXIT_FAILURE);
01730 }
01731
01732 strncpy( chLocal, chInputLine, INPUT_LINE_LENGTH );
01733
01734
01735 if( chLocal[INPUT_LINE_LENGTH-1] != '\0' )
01736 {
01737 chLocal[INPUT_LINE_LENGTH-1] = '\0';
01738 fprintf(ioQQQ," PROBLEM cdRead, while parsing the following input line:\n %s\n",
01739 chInputLine);
01740 fprintf(ioQQQ," found that the line is longer than %i characters, the longest possible line.\n",
01741 INPUT_LINE_LENGTH-1);
01742 fprintf(ioQQQ," Please make the command line no longer than this limit.\n");
01743 }
01744
01745
01746
01747 if( (chEOL = strchr_s(chLocal, '\n' ) ) != NULL )
01748 {
01749 *chEOL = '\0';
01750 }
01751 if( (chEOL = strchr_s(chLocal, '%' ) ) != NULL )
01752 {
01753 *chEOL = '\0';
01754 }
01755
01756 if( (chEOL = strchr_s(chLocal, '#' ) ) != NULL )
01757 {
01758 *chEOL = '\0';
01759 }
01760 if( (chEOL = strchr_s(chLocal, ';' ) ) != NULL )
01761 {
01762 *chEOL = '\0';
01763 }
01764 if( (chEOL = strstr_s(chLocal, "//" ) ) != NULL )
01765 {
01766 *chEOL = '\0';
01767 }
01768
01769
01770
01771 if( (chEOL = strchr_s( chLocal, '\0' )) == NULL )
01772 TotalInsanity();
01773
01774
01775
01776 if( chEOL-chLocal < INPUT_LINE_LENGTH-2 )
01777 {
01778 strcat( chLocal, " " );
01779 }
01780
01781
01782 strcpy( input.chCardSav[input.nSave], chLocal );
01783
01784
01785 strcpy( chCARD, chLocal );
01786 caps(chCARD);
01787
01788
01789 lgComment = false;
01790 if( strncmp(chCARD , "C ", 2 )==0 )
01791 lgComment = true;
01792
01793
01794 if( strncmp(chCARD,"TRACE",5) == 0 )
01795 trace.lgTrace = true;
01796
01797
01798 if( trace.lgTrace )
01799 fprintf( ioQQQ,"cdRead=%s=\n",input.chCardSav[input.nSave] );
01800
01801 {
01802
01803 char chFilename[INPUT_LINE_LENGTH],
01804 chTemp[INPUT_LINE_LENGTH];
01805
01806 strcpy( chTemp , input.chCardSav[input.nSave] );
01807 GetQuote( chFilename , chCARD , chTemp , false );
01808 }
01809
01810
01811 if( !lgComment && nMatch("VARY",chCARD) )
01812
01813 optimize.lgVaryOn = true;
01814
01815
01816 if( strncmp(chCARD,"NO VARY",7) == 0 )
01817 optimize.lgNoVary = true;
01818
01819
01820 if( strncmp(chCARD,"GRID",4) == 0 )
01821 {
01822 grid.lgGrid = true;
01823 ++grid.nGridCommands;
01824 }
01825
01826
01827
01828 if( strncmp(chCARD,"NO BUFF",7) == 0 )
01829 {
01830
01831
01832
01833 FILE *test = stdout;
01834 if( ioQQQ == test )
01835 {
01836 fprintf( ioQQQ, " cdRead found NO BUFFERING command, redirecting output to stderr now.\n" );
01837
01838 fflush( ioQQQ );
01839
01840 ioQQQ = stderr;
01841
01842 input.lgSetNoBuffering = false;
01843 }
01844 }
01845
01846
01847 if( strncmp(chCARD,"OPTI",4) == 0 || strncmp(chCARD,"GRID",4) == 0 )
01848
01849 optimize.lgOptimr = true;
01850
01851 return NKRD - input.nSave;
01852 }
01853
01854
01855 void cdClosePunchFiles()
01856 {
01857
01858 DEBUG_ENTRY( "cdClosePunchFiles()" );
01859
01860 CloseSaveFiles( true );
01861 return;
01862 }