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