00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "trace.h"
00009 #include "hydro_bauman.h"
00010 #include "hydroeinsta.h"
00011 #include "thirdparty.h"
00012 #include "phycon.h"
00013 #include "atmdat.h"
00014 #include "taulines.h"
00015 #include "thermal.h"
00016 #include "thirdparty.h"
00017 #include "abund.h"
00018 #include "rt.h"
00019 #include "continuum.h"
00020 #include "parser.h"
00021 #include "physconst.h"
00022
00023
00024 STATIC void dgaunt(void);
00025
00026
00027 STATIC void DrvHyas(void);
00028
00029
00030 STATIC void DrvEscP( void );
00031
00032
00033 STATIC void DrvCaseBHS(void);
00034
00035 void ParseDrive(Parser &p )
00036 {
00037 bool lgEOL;
00038 long int n,
00039 i;
00040 double fac,
00041 zed;
00042
00043 DEBUG_ENTRY( "ParseDrive()" );
00044
00045
00046
00047
00048 if( p.nMatch("FFMT") )
00049 {
00050
00051 char chInput[INPUT_LINE_LENGTH];
00052 fprintf( ioQQQ, " FFmtRead ParseDrive entered. Enter number.\n" );
00053 lgEOL = false;
00054 while( !lgEOL )
00055 {
00056 if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
00057 {
00058 fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
00059 cdEXIT(EXIT_FAILURE);
00060 }
00061 i = 1;
00062 fac = FFmtRead(chInput,&i,sizeof(chInput),&lgEOL);
00063 if( lgEOL )
00064 {
00065 fprintf( ioQQQ, " FFmtRead hit the EOL with no value, return=%10.2e\n",
00066 fac );
00067 break;
00068 }
00069 else if( fac == 0. )
00070 {
00071 break;
00072 }
00073 else
00074 {
00075 fprintf( ioQQQ, " FFmtRead returned with value%11.4e\n",
00076 fac );
00077 }
00078 fprintf( ioQQQ, " Enter 0 to stop, or another value.\n" );
00079 }
00080 fprintf( ioQQQ, " FFmtRead ParseDrive exits.\n" );
00081 }
00082
00083 else if( p.nMatch("CASE") )
00084 {
00085
00086 DrvCaseBHS( );
00087 }
00088
00089 else if( p.nMatch("CDLI") )
00090 {
00091
00092 trace.lgDrv_cdLine = true;
00093 }
00094
00095 else if( p.nMatch(" E1 ") )
00096 {
00097
00098
00099 double tau = 1.275;
00100 for( i=0; i<50; ++i )
00101 {
00102 fprintf(ioQQQ,"tau\t%.3e\t exp-tau\t%.5e\t e1 tau\t%.5e \t e2 "
00103 "\t%.5e \te2n %.5e \t e3\t%.5e \t e4\t%.5e \n",
00104 tau , sexp(tau) , ee1(tau) , e2(tau ), expn(2, tau) ,
00105 expn(3 , tau ), expn(4 , tau ) );
00106 tau = pow( 10. , ((double)i/4. - 9.) );
00107 }
00108 cdEXIT(EXIT_SUCCESS);
00109 }
00110
00111 else if( p.nMatch("ESCA") )
00112 {
00113
00114 DrvEscP( );
00115 }
00116
00117 else if( p.nMatch("HYAS") )
00118 {
00119
00120 DrvHyas();
00121 }
00122
00123 else if( p.nMatch("GAUN") )
00124 {
00125
00126 dgaunt();
00127 }
00128
00129 else if( p.nMatch("POIN") )
00130 {
00131
00132 fprintf( ioQQQ, " Define entire model first, later I will ask for energies.\n" );
00133 trace.lgPtrace = true;
00134 }
00135
00136 else if( p.nMatch("PUMP") )
00137 {
00138 char chInput[INPUT_LINE_LENGTH];
00139 lgEOL = false;
00140 fprintf( ioQQQ, " Continuum pump ParseDrive entered - Enter log tau\n" );
00141 while( !lgEOL )
00142 {
00143 if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
00144 {
00145 fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
00146 cdEXIT(EXIT_FAILURE);
00147 }
00148
00149 i = 1;
00150 fac = FFmtRead(chInput,&i,sizeof(chInput),&lgEOL);
00151 if( lgEOL )
00152 break;
00153 fac = pow(10.,fac);
00154 fprintf( ioQQQ, " Tau =%11.4e\n", fac );
00155 (*TauDummy).Emis().TauIn() = (realnum)fac;
00156
00157
00158 fac = DrvContPump(*TauDummy,1.f);
00159 fprintf( ioQQQ, " ContPump =%11.4e\n", fac );
00160 fprintf( ioQQQ, " Enter null to stop, or another value.\n" );
00161 }
00162 fprintf( ioQQQ, " ContPump ParseDrive exits.\n" );
00163 }
00164
00165 else if( p.nMatch("STAR") )
00166 {
00167 char chInput[INPUT_LINE_LENGTH];
00168
00169 for( i=0; i < 40; i++ )
00170 {
00171 zed = ((double)i+1.)/4. + 0.01;
00172 sprintf( chInput, "starburst, zed=%10.4f", zed );
00173 p.setline(chInput);
00174 abund_starburst(p);
00175 fprintf( ioQQQ, "%8.1e", zed );
00176 for(n=0; n < LIMELM; n++)
00177 fprintf( ioQQQ, "%8.1e", abund.solar[n] );
00178 fprintf( ioQQQ, "\n" );
00179 }
00180 }
00181
00182 else if( p.nMatch("VOIGT") )
00183 {
00184
00185 FILE *ioVOIGT = fopen("voigt.dat" , "w");
00186 fprintf(ioVOIGT,"x \\ a");
00187 const realnum DampLogMin = -4., DampLogMax = 4.01;
00188 for( realnum damplog=DampLogMin; damplog<DampLogMax; ++damplog)
00189 fprintf(ioVOIGT,"\tlog a=%.3e",pow(10.,damplog));
00190 fprintf(ioVOIGT , "\n");
00191
00192 for( realnum x=-2.; x<5.;x+=0.05)
00193 {
00194 realnum xlin = (realnum)pow(10.,x);
00195 fprintf(ioVOIGT,"%.3e",xlin);
00196 for( realnum damplog=DampLogMin; damplog<DampLogMax; ++damplog)
00197 {
00198 realnum xval[1];
00199 xval[0] = xlin;
00200 realnum damp = (realnum)pow(10. , damplog);
00201 realnum yval[1];
00202 VoigtH(damp,xval,yval,1);
00203 fprintf(ioVOIGT , "\t%.3e",yval[0]);
00204 }
00205 fprintf(ioVOIGT , "\n");
00206 }
00207 fclose(ioVOIGT);
00208 cdEXIT(EXIT_SUCCESS);
00209 }
00210
00211 else
00212 {
00213 fprintf( ioQQQ,
00214 " Unrecognized key; keys are CASE, CDLIne, E1 , ESCApe, FFMTread, GAUNt, "
00215 "HYAS, POINt, PUMP, STAR, and VOIGt. Sorry.\n" );
00216 cdEXIT(EXIT_FAILURE);
00217 }
00218 return;
00219 }
00220
00221
00222 STATIC void DrvEscP( void )
00223 {
00224 char chCard[INPUT_LINE_LENGTH];
00225 bool lgEOL;
00226 long i;
00227 double tau;
00228
00229 DEBUG_ENTRY( "DrvEscP()" );
00230
00231
00232
00233 fprintf( ioQQQ, " Enter the log of the one-sided optical depth; line with no number to stop.\n" );
00234
00235 lgEOL = false;
00236 while( !lgEOL )
00237 {
00238 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00239 {
00240 break;
00241 }
00242
00243 i = 1;
00244 tau = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00245 if( lgEOL )
00246 {
00247 break;
00248 }
00249
00250 tau = pow(10.,tau);
00251 fprintf( ioQQQ, "tau was %e\n", tau );
00252 fprintf( ioQQQ, " ESCINC=%11.3e\n", esc_PRD_1side(tau,1e-4) );
00253 fprintf( ioQQQ, " ESCCOM=%11.3e\n", esc_CRDwing_1side(tau,1e-4 ) );
00254 fprintf( ioQQQ, " ESCA0K2=%11.3e\n", esca0k2(tau) );
00255
00256 }
00257 return;
00258 }
00259
00260
00261 STATIC void DrvCaseBHS(void)
00262 {
00263 char chCard[INPUT_LINE_LENGTH];
00264 bool lgEOL,
00265 lgHit;
00266 long int i,
00267 n1,
00268 nelem ,
00269 n2;
00270 double Temperature,
00271 Density;
00272
00273 DEBUG_ENTRY( "DrvCaseBHS()" );
00274
00275
00276
00277
00278
00279 fprintf(ioQQQ," I will get needed H data files. This will take a second.\n");
00280 atmdat_readin();
00281
00282 {
00283
00284
00285 enum {DEBUG_LOC=false};
00286
00287 if( DEBUG_LOC )
00288 {
00289 double xLyman , alpha;
00290 long int ipHi;
00291 nelem = 2;
00292 Temperature = 2e4;
00293 Density = 1e2;
00294 for( ipHi=3; ipHi<25; ++ipHi )
00295 {
00296 double photons = (1./POW2(ipHi-1.)-1./POW2((double)ipHi) ) /(1.-1./ipHi/ipHi );
00297 xLyman = atmdat_HS_caseB(1,ipHi, nelem,Temperature , Density , 'A' );
00298 alpha = atmdat_HS_caseB(ipHi-1,ipHi, nelem,Temperature , Density , 'A' );
00299 fprintf(ioQQQ,"%li\t%.3e\t%.3e\n", ipHi, xLyman/alpha*photons, photons );
00300 }
00301 cdEXIT(EXIT_SUCCESS);
00302 }
00303 }
00304
00305
00306 lgHit = false;
00307 nelem = 0;
00308 while( !lgHit )
00309 {
00310 fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" );
00311 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00312 {
00313 fprintf( ioQQQ, " error getting species \n" );
00314 cdEXIT(EXIT_FAILURE);
00315 }
00316
00317 i = 1;
00318 nelem = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00319 if( lgEOL || nelem< 1 || nelem > 2 )
00320 {
00321 fprintf( ioQQQ, " must be either 1 or 2!\n" );
00322 }
00323 else
00324 {
00325 lgHit = true;
00326 }
00327 }
00328
00329 fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n");
00330 fprintf(ioQQQ," The density is always a log.\n");
00331 fprintf(ioQQQ," The order of the quantum numbers do not matter.\n");
00332 fprintf(ioQQQ," The smallest must not be smaller than 2,\n");
00333 fprintf(ioQQQ," and the largest must not be larger than 25.\n");
00334 fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n");
00335 fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n");
00336
00337 lgHit = true;
00338
00339 while( lgHit )
00340 {
00341 fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" );
00342 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00343 {
00344 fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
00345 cdEXIT(EXIT_FAILURE);
00346 }
00347
00348 i = 1;
00349 Temperature = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00350 if( lgEOL )
00351 {
00352 fprintf( ioQQQ, " error getting temperature!\n" );
00353 break;
00354 }
00355
00356
00357 if( Temperature < 10. )
00358 {
00359 Temperature = pow(10., Temperature );
00360 }
00361 fprintf(ioQQQ," Temperature is %g\n", Temperature );
00362
00363 Density = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00364 if( lgEOL )
00365 {
00366 fprintf( ioQQQ, " error getting density!\n" );
00367 break;
00368 }
00369 Density = pow(10., Density );
00370 fprintf(ioQQQ," Density is %g\n", Density );
00371
00372
00373 n1 = (long)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00374 if( lgEOL )
00375 {
00376 fprintf( ioQQQ, " error getting quantum number!\n" );
00377 break;
00378 }
00379
00380 n2 = (long)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00381 if( lgEOL )
00382 {
00383 fprintf( ioQQQ, " error getting quantum number!\n" );
00384 break;
00385 }
00386
00387 if( MAX2( n1 , n2 ) > 25 )
00388 {
00389 fprintf( ioQQQ," The limits of the HS tables are 2 <= n <= 25. Sorry.\n");
00390 break;
00391 }
00392
00393 fprintf( ioQQQ,
00394 " 4pJ(%ld,%ld)/n_e n_p=%11.3e\n",
00395 n1, n2,
00396 atmdat_HS_caseB(n1,n2, nelem,Temperature , Density , 'B' ) );
00397
00398
00399 #if 0
00400 {
00401 long j;
00402 double tempTable[33] = {
00403 11870.,12490.,12820.,
00404 11060.,17740.,12560.,
00405 16390.,16700.,11360.,
00406 10240.,20740.,12030.,
00407 14450.,19510.,12550.,
00408 16470.,16560.,12220.,
00409 15820.,12960.,10190.,
00410 12960.,14060.,12560.,
00411 11030.,10770.,13360.,
00412 10780.,11410.,11690.,
00413 12500.,13190.,21120. };
00414 double edenTable[33] = {
00415 10.,270.,80.,10.,70.,
00416 110.,200.,10.,40.,90.,
00417 340.,80.,60.,340.,30.,
00418 120.,10.,50.,450.,30.,
00419 180.,20.,170.,60.,20.,
00420 40.,30.,20.,100.,130.,
00421 10.,10.,110. };
00422
00423
00424 for( j=0; j<33; j++ )
00425 {
00426 double halpha, hbeta, hgamma;
00427
00428 halpha = atmdat_HS_caseB(2,3, 1,tempTable[j] , edenTable[j] , 'B' );
00429 hbeta = atmdat_HS_caseB(2,4, 1,tempTable[j] , edenTable[j] , 'B' );
00430 hgamma = atmdat_HS_caseB(2,5, 1,tempTable[j] , edenTable[j] , 'B' );
00431
00432 fprintf( ioQQQ, "%e\t%e\t%e\t%e\n",
00433 tempTable[j],
00434 edenTable[j],
00435 halpha/hbeta,
00436 hgamma/hbeta );
00437 }
00438 }
00439 #endif
00440 }
00441
00442 fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
00443 cdEXIT(EXIT_FAILURE);
00444
00445 }
00446
00447
00448 STATIC void DrvHyas(void)
00449 {
00450 char chCard[INPUT_LINE_LENGTH];
00451 bool lgEOL;
00452 long int i, nHi, lHi, nLo, lLo;
00453
00454 DEBUG_ENTRY( "DrvHyas()" );
00455
00456
00457
00458
00459 nHi = 1;
00460
00461 while( nHi != 0 )
00462 {
00463 fprintf( ioQQQ, " Enter four quantum numbers (n, l, n', l'), null line to stop.\n" );
00464 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00465 {
00466 fprintf( ioQQQ, " error getting drvhyas \n" );
00467 cdEXIT(EXIT_FAILURE);
00468 }
00469
00470 i = 1;
00471 nHi = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00472 if( lgEOL )
00473 break;
00474
00475 lHi = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00476 if( lgEOL )
00477 {
00478 fprintf( ioQQQ, " must be four numbers!\n" );
00479 break;
00480 }
00481
00482 if( lHi >= nHi )
00483 {
00484 fprintf( ioQQQ, " l must be less than n!\n" );
00485 break;
00486 }
00487
00488 nLo = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00489 if( lgEOL )
00490 {
00491 fprintf( ioQQQ, " must be four numbers!\n" );
00492 break;
00493 }
00494
00495 lLo = (long int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00496 if( lgEOL )
00497 {
00498 fprintf( ioQQQ, " must be four numbers!\n" );
00499 break;
00500 }
00501
00502 if( lLo >= nLo )
00503 {
00504 fprintf( ioQQQ, " l must be less than n!\n" );
00505 break;
00506 }
00507
00508 if( nLo > nHi )
00509 {
00510 long nTemp, lTemp;
00511
00512
00513 nTemp = nLo;
00514 lTemp = lLo;
00515 nLo = nHi;
00516 lLo = lHi;
00517 nHi = nTemp;
00518 lHi = lTemp;
00519 }
00520
00521 fprintf( ioQQQ, " A(%3ld,%3ld->%3ld,%3ld)=%11.3e\n",
00522 nHi, lHi, nLo, lLo,
00523 H_Einstein_A( nHi, lHi, nLo, lLo, 1 ) );
00524
00525 }
00526 fprintf( ioQQQ, " Driver exits, enter next line.\n" );
00527
00528 return;
00529 }
00530
00531
00532 STATIC void dgaunt(void)
00533 {
00534 char chCard[INPUT_LINE_LENGTH];
00535 bool lgEOL;
00536 int inputflag;
00537 long int i,
00538 ierror;
00539 realnum enerlin[1];
00540 double SaveTemp;
00541 double z,mygaunt=0.;
00542 double loggamma2, logu;
00543
00544 DEBUG_ENTRY( "dgaunt()" );
00545
00546 SaveTemp = phycon.te;
00547
00548
00549
00550
00551 fprintf( ioQQQ, " Enter 0 to input temp, energy, and net charge, or 1 for u and gamma**2.\n" );
00552 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00553 {
00554 fprintf( ioQQQ, " dgaunt error getting magic number\n" );
00555 cdEXIT(EXIT_FAILURE);
00556 }
00557 i = 1;
00558 inputflag = (int)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00559
00560 if( inputflag == 0 )
00561 {
00562 fprintf( ioQQQ, " Enter the temperature (log if <=10), energy (Ryd), and net charge. Null line to stop.\n" );
00563
00564
00565 ierror = 0;
00566 while( ierror == 0 )
00567 {
00568 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00569 {
00570 fprintf( ioQQQ, " dgaunt error getting magic number\n" );
00571 cdEXIT(EXIT_FAILURE);
00572 }
00573 i = 1;
00574 phycon.alogte = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00575
00576 if( lgEOL )
00577 {
00578 fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
00579 break;
00580 }
00581
00582 double TeNew;
00583 if( phycon.alogte > 10. )
00584 {
00585 TeNew = phycon.alogte;
00586 }
00587 else
00588 {
00589 TeNew = pow(10.,phycon.alogte);
00590 }
00591 TempChange(TeNew , false);
00592
00593 enerlin[0] = (realnum)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00594 if( lgEOL || enerlin[0] == 0. )
00595 {
00596 fprintf( ioQQQ, " Sorry, but there should be two more numbers, energy and charge.\n" );
00597 }
00598
00599 z = (double)FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00600 if( lgEOL || z == 0. )
00601 {
00602 fprintf( ioQQQ, " Sorry, but there should be a third number, charge.\n" );
00603 }
00604
00605
00606 mygaunt = cont_gaunt_calc( (double)phycon.te, z, enerlin[0] );
00607
00608 fprintf( ioQQQ, " Using my routine, Gff= \t" );
00609 fprintf( ioQQQ, "%11.3e\n", mygaunt );
00610
00611 }
00612 }
00613 else
00614 {
00615
00616
00617
00618 fprintf( ioQQQ, " Enter log u and log gamma2. Null line to stop.\n" );
00619 ierror = 0;
00620 while( ierror == 0 )
00621 {
00622 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00623 {
00624 fprintf( ioQQQ, " dgaunt error getting magic number\n" );
00625 cdEXIT(EXIT_FAILURE);
00626 }
00627 i = 1;
00628 logu = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00629
00630 if( lgEOL )
00631 {
00632 fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
00633 break;
00634 }
00635
00636 loggamma2 = FFmtRead(chCard,&i,sizeof(chCard),&lgEOL);
00637 if( lgEOL )
00638 {
00639 fprintf( ioQQQ, " Sorry, but there should be another numbers, log gamma2.\n" );
00640 }
00641
00642
00643 mygaunt = cont_gaunt_calc( TE1RYD/pow(10.,loggamma2), 1., pow(10.,logu-loggamma2) );
00644
00645 TempChange(TE1RYD/pow(10.,loggamma2) , false);
00646
00647 fprintf( ioQQQ, " Using my routine, Gaunt factor is" );
00648 fprintf( ioQQQ, "%11.3e\n", mygaunt );
00649 }
00650 }
00651
00652 TempChange(SaveTemp , false);
00653 return;
00654 }