00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "trace.h"
00009 #include "hydroeinsta.h"
00010 #include "phycon.h"
00011 #include "atmdat.h"
00012 #include "taulines.h"
00013 #include "thermal.h"
00014 #include "abund.h"
00015 #include "rt.h"
00016 #include "continuum.h"
00017 #include "parse.h"
00018 #include "physconst.h"
00019
00020
00021 static void dgaunt(void);
00022
00023
00024 static void DrvHyas(void);
00025
00026
00027 static void DrvEscP( void );
00028
00029
00030 static void DrvCaseBHS(void);
00031
00032 void ParseDrive(char *chCard )
00033 {
00034 bool lgEOL;
00035 long int n,
00036 i;
00037 double fac,
00038 zed;
00039
00040 DEBUG_ENTRY( "ParseDrive()" );
00041
00042
00043
00044
00045 if( nMatch("FFMT",chCard) )
00046 {
00047
00048 char chInput[INPUT_LINE_LENGTH];
00049 fprintf( ioQQQ, " FFmtRead ParseDrive entered. Enter number.\n" );
00050 lgEOL = false;
00051 while( !lgEOL )
00052 {
00053 if( fgets( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
00054 {
00055 fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
00056 puts( "[Stop in ParseDrive]" );
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(1);
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( fgets( chInput , (int)sizeof(chInput) , ioStdin ) == NULL )
00138 {
00139 fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" );
00140 puts( "[Stop in ParseDrive]" );
00141 cdEXIT(EXIT_FAILURE);
00142 }
00143
00144 i = 1;
00145 fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL);
00146 if( lgEOL )
00147 break;
00148 fac = pow(10.,fac);
00149 fprintf( ioQQQ, " Tau =%11.4e\n", fac );
00150 TauDummy.TauIn = (float)fac;
00151 fac = DrvContPump(&TauDummy);
00152 fprintf( ioQQQ, " ContPump =%11.4e\n", fac );
00153 fprintf( ioQQQ, " Enter null to stop, or another value.\n" );
00154 }
00155 fprintf( ioQQQ, " ContPump ParseDrive exits.\n" );
00156 }
00157
00158 else if( nMatch("STAR",chCard) )
00159 {
00160 char chInput[INPUT_LINE_LENGTH];
00161
00162 for( i=0; i < 40; i++ )
00163 {
00164 zed = (float)(i+1)/4. + 0.01;
00165 sprintf( chInput, "starburst, zed=%10.4f", zed );
00166 abund_starburst(chInput);
00167 fprintf( ioQQQ, "%8.1e", zed );
00168 for(n=0; n < LIMELM; n++)
00169 fprintf( ioQQQ, "%8.1e", abund.solar[n] );
00170 fprintf( ioQQQ, "\n" );
00171 }
00172 }
00173
00174 else
00175 {
00176 fprintf( ioQQQ,
00177 " Unrecognized key; keys are FFmtRead, CASE, GAUNT, HYAS, POINT, MOLE, STAR, "
00178 "PUMP, and ESCApe. Sorry.\n" );
00179 puts( "[Stop in ParseDrive]" );
00180 cdEXIT(EXIT_FAILURE);
00181 }
00182
00183 DEBUG_EXIT( "ParseDrive()" );
00184 return;
00185 }
00186
00187
00188 static void DrvEscP( void )
00189 {
00190 char chCard[INPUT_LINE_LENGTH];
00191 bool lgEOL;
00192 long i;
00193 double tau;
00194
00195 DEBUG_ENTRY( "DrvEscP()" );
00196
00197
00198
00199 fprintf( ioQQQ, " Enter the log of the one-sided optical depth; line with no number to stop.\n" );
00200
00201 lgEOL = false;
00202 while( !lgEOL )
00203 {
00204 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00205 {
00206 break;
00207 }
00208
00209 i = 1;
00210 tau = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00211 if( lgEOL )
00212 {
00213 break;
00214 }
00215
00216 tau = pow(10.,tau);
00217 fprintf( ioQQQ, "tau was %e\n", tau );
00218 fprintf( ioQQQ, " ESCINC=%11.3e\n", esc_PRD_1side(tau,1e-4) );
00219 fprintf( ioQQQ, " ESCCOM=%11.3e\n", esc_CRDwing_1side(tau,1e-4 ) );
00220 fprintf( ioQQQ, " ESCA0K2=%11.3e\n", esca0k2(tau) );
00221
00222 }
00223
00224 DEBUG_EXIT( "DrvEscP()" );
00225 return;
00226 }
00227
00228
00229 static void DrvCaseBHS(void)
00230 {
00231 char chCard[INPUT_LINE_LENGTH];
00232 bool lgEOL,
00233 lgHit;
00234 long int i,
00235 n1,
00236 nelem ,
00237 n2;
00238 double Temperature,
00239 Density;
00240
00241 DEBUG_ENTRY( "DrvCaseBHS()" );
00242
00243
00244
00245
00246
00247 fprintf(ioQQQ," I will get needed H data files. This will take a second.\n");
00248 atmdat_readin();
00249
00250 {
00251
00252
00253 enum {DEBUG_LOC=false};
00254
00255 if( DEBUG_LOC )
00256 {
00257 double xLyman , alpha;
00258 long int ipHi;
00259 nelem = 2;
00260 Temperature = 2e4;
00261 Density = 1e2;
00262 for( ipHi=3; ipHi<25; ++ipHi )
00263 {
00264 double photons = (1./POW2(ipHi-1.)-1./POW2((double)ipHi) ) /(1.-1./ipHi/ipHi );
00265 xLyman = atmdat_HS_caseB(1,ipHi, nelem,Temperature , Density , 'A' );
00266 alpha = atmdat_HS_caseB(ipHi-1,ipHi, nelem,Temperature , Density , 'A' );
00267 fprintf(ioQQQ,"%li\t%.3e\t%.3e\n", ipHi, xLyman/alpha*photons, photons );
00268 }
00269 cdEXIT(EXIT_SUCCESS);
00270 }
00271 }
00272
00273
00274 lgHit = false;
00275 nelem = 0;
00276 while( !lgHit )
00277 {
00278 fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" );
00279 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00280 {
00281 fprintf( ioQQQ, " error getting species \n" );
00282 puts( "[Stop in DrvCaseBHS]" );
00283 cdEXIT(EXIT_FAILURE);
00284 }
00285
00286 i = 1;
00287 nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00288 if( lgEOL || nelem< 1 || nelem > 2 )
00289 {
00290 fprintf( ioQQQ, " must be either 1 or 2!\n" );
00291 }
00292 else
00293 {
00294 lgHit = true;
00295 }
00296 }
00297
00298 fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n");
00299 fprintf(ioQQQ," The density is always a log.\n");
00300 fprintf(ioQQQ," The order of the quantum numbers do not matter.\n");
00301 fprintf(ioQQQ," The smallest must not be smaller than 2,\n");
00302 fprintf(ioQQQ," and the largest must not be larger than 25.\n");
00303 fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n");
00304 fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n");
00305
00306 lgHit = true;
00307
00308 while( lgHit )
00309 {
00310 fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" );
00311 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00312 {
00313 fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" );
00314 puts( "[Stop in DrvCaseBHS]" );
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 puts( "[Stop in DrvCaseBHS]" );
00414 cdEXIT(EXIT_FAILURE);
00415
00416 }
00417
00418
00419 static void DrvHyas(void)
00420 {
00421 char chCard[INPUT_LINE_LENGTH];
00422 bool lgEOL;
00423 long int i,
00424 ihi,
00425 low;
00426
00427 DEBUG_ENTRY( "DrvHyas()" );
00428
00429
00430
00431
00432 ihi = 1;
00433
00434 while( ihi != 0 )
00435 {
00436 fprintf( ioQQQ, " Enter two quantum numbers, null line to stop.\n" );
00437 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00438 {
00439 fprintf( ioQQQ, " error getting drvhyas \n" );
00440 puts( "[Stop in drvhyas]" );
00441 cdEXIT(EXIT_FAILURE);
00442 }
00443
00444 i = 1;
00445 low = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00446 if( lgEOL )
00447 break;
00448
00449 ihi = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00450 if( lgEOL )
00451 {
00452 fprintf( ioQQQ, " must be two numbers!\n" );
00453 break;
00454 }
00455 fprintf( ioQQQ, " A(%3ld,%3ld)=%11.3e\n",
00456 MIN2(low,ihi), MAX2(low,ihi), HydroEinstA(low,ihi) );
00457
00458 }
00459 fprintf( ioQQQ, " Driver exits, enter next line.\n" );
00460
00461
00462 DEBUG_EXIT( "DrvHyas()" );
00463 return;
00464 }
00465
00466
00467 static void dgaunt(void)
00468 {
00469 char chCard[INPUT_LINE_LENGTH];
00470 bool lgEOL;
00471 int inputflag;
00472 long int i,
00473 ierror;
00474 float enerlin[1];
00475 double SaveTemp;
00476 double z,mygaunt=0.;
00477 double loggamma2, logu;
00478
00479 DEBUG_ENTRY( "dgaunt()" );
00480
00481 SaveTemp = phycon.te;
00482
00483
00484
00485
00486 fprintf( ioQQQ, " Enter 0 to input temp, energy, and net charge, or 1 for u and gamma**2.\n" );
00487 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00488 {
00489 fprintf( ioQQQ, " dgaunt error getting magic number\n" );
00490 puts( "[Stop in dgaunt]" );
00491 cdEXIT(EXIT_FAILURE);
00492 }
00493 i = 1;
00494 inputflag = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00495
00496 if( inputflag == 0 )
00497 {
00498 fprintf( ioQQQ, " Enter the temperature (log if <=10), energy (Ryd), and net charge. Null line to stop.\n" );
00499
00500
00501 ierror = 0;
00502 while( ierror == 0 )
00503 {
00504 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00505 {
00506 fprintf( ioQQQ, " dgaunt error getting magic number\n" );
00507 puts( "[Stop in dgaunt]" );
00508 cdEXIT(EXIT_FAILURE);
00509 }
00510 i = 1;
00511 phycon.alogte = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00512
00513 if( lgEOL )
00514 {
00515 fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
00516 break;
00517 }
00518
00519 if( phycon.alogte > 10. )
00520 {
00521 phycon.te = phycon.alogte;
00522 phycon.alogte = log10(phycon.te);
00523 }
00524 else
00525 {
00526 phycon.te = pow(10.,phycon.alogte);
00527 }
00528 tfidle(false);
00529
00530 enerlin[0] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00531 if( lgEOL || enerlin[0] == 0. )
00532 {
00533 fprintf( ioQQQ, " Sorry, but there should be two more numbers, energy and charge.\n" );
00534 }
00535
00536 z = (double)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00537 if( lgEOL || z == 0. )
00538 {
00539 fprintf( ioQQQ, " Sorry, but there should be a third number, charge.\n" );
00540 }
00541
00542
00543 mygaunt = cont_gaunt_calc( (double)phycon.te, z, enerlin[0] );
00544
00545 fprintf( ioQQQ, " Using my routine, Gff= \t" );
00546 fprintf( ioQQQ, "%11.3e\n", mygaunt );
00547
00548 }
00549 }
00550 else
00551 {
00552
00553
00554
00555 fprintf( ioQQQ, " Enter log u and log gamma2. Null line to stop.\n" );
00556 ierror = 0;
00557 while( ierror == 0 )
00558 {
00559 if( fgets( chCard , (int)sizeof(chCard) , ioStdin ) == NULL )
00560 {
00561 fprintf( ioQQQ, " dgaunt error getting magic number\n" );
00562 puts( "[Stop in dgaunt]" );
00563 cdEXIT(EXIT_FAILURE);
00564 }
00565 i = 1;
00566 logu = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00567
00568 if( lgEOL )
00569 {
00570 fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" );
00571 break;
00572 }
00573
00574 loggamma2 = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00575 if( lgEOL )
00576 {
00577 fprintf( ioQQQ, " Sorry, but there should be another numbers, log gamma2.\n" );
00578 }
00579
00580
00581 mygaunt = cont_gaunt_calc( TE1RYD/pow(10.,loggamma2), 1., pow(10.,logu-loggamma2) );
00582
00583 phycon.te = (TE1RYD/pow(10.,loggamma2));
00584 tfidle(false);
00585
00586 fprintf( ioQQQ, " Using my routine, Gaunt factor is" );
00587 fprintf( ioQQQ, "%11.3e\n", mygaunt );
00588 }
00589 }
00590
00591 phycon.te = SaveTemp;
00592 tfidle(false);
00593
00594 DEBUG_EXIT( "dgaunt()" );
00595 return;
00596 }
00597