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