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