00001
00002
00003
00004
00005
00006 #include "cddefines.h"
00007 #include "physconst.h"
00008 #include "opacity.h"
00009 #include "iso.h"
00010 #include "dense.h"
00011 #include "phycon.h"
00012 #include "stopcalc.h"
00013 #include "continuum.h"
00014 #include "trace.h"
00015 #include "rfield.h"
00016 #include "doppvel.h"
00017 #include "radius.h"
00018 #include "wind.h"
00019 #include "thermal.h"
00020
00021
00022 STATIC void tauff(void);
00023
00024 STATIC void velset(void);
00025
00026 STATIC void FillGFF(void);
00027
00028 STATIC realnum InterpolateGff( long charge , double ERyd );
00029 STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx);
00030 STATIC int LinterpVector(realnum **t, realnum *v, long lta , long ltb, realnum *yy , long ny, realnum **a);
00031 STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j);
00032
00036 STATIC void tfidle(
00037 bool lgForceUpdate);
00038
00039 static long lgGffNotFilled = true;
00040
00041 const long N_TE_GFF = 21;
00042 const long N_PHOTON_GFF = 145;
00043 static realnum ***GauntFF;
00044 static realnum **GauntFF_T;
00045
00046 static realnum TeGFF[N_TE_GFF];
00047
00048 static realnum PhoGFF[N_PHOTON_GFF];
00049
00052 void TempChange(
00053 double TempNew ,
00054
00055 bool lgForceUpdate)
00056 {
00057
00058 DEBUG_ENTRY( "TempChange()" );
00059
00060
00061 if( TempNew >= phycon.TEMP_LIMIT_HIGH )
00062 {
00063
00064 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
00065 " is above the upper limit of the code, %.3eK.\n",
00066 TempNew , phycon.TEMP_LIMIT_HIGH );
00067 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
00068
00069 TempNew = phycon.TEMP_LIMIT_HIGH*0.999;
00070 lgAbort = true;
00071 }
00072 else if( TempNew <= phycon.TEMP_LIMIT_LOW )
00073 {
00074
00075 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
00076 " is below the lower limit of the code, %.3eK.\n",
00077 TempNew , phycon.TEMP_LIMIT_LOW );
00078 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
00079 TempNew = phycon.TEMP_LIMIT_LOW*1.0001;
00080 lgAbort = true;
00081 }
00082 else if( TempNew < StopCalc.TeFloor )
00083 {
00084
00085
00086
00087 thermal.lgTemperatureConstant = true;
00088 thermal.ConstTemp = (realnum)StopCalc.TeFloor;
00089 phycon.te = thermal.ConstTemp;
00090
00091
00092 }
00093 else
00094 {
00095
00096 phycon.te = TempNew;
00097 }
00098
00099
00100 tfidle( lgForceUpdate );
00101 return;
00102 }
00106 void TempChange(
00107 double TempNew )
00108 {
00109
00110 DEBUG_ENTRY( "TempChange()" );
00111
00112
00113 if( TempNew >= phycon.TEMP_LIMIT_HIGH )
00114 {
00115
00116 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
00117 " is above the upper limit of the code, %.3eK.\n",
00118 TempNew , phycon.TEMP_LIMIT_HIGH );
00119 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
00120
00121 TempNew = phycon.TEMP_LIMIT_HIGH*0.999;
00122 }
00123 else if( TempNew <= phycon.TEMP_LIMIT_LOW )
00124 {
00125
00126 fprintf(ioQQQ," PROBLEM DISASTER - the kinetic temperature, %.3eK,"
00127 " is below the lower limit of the code, %.3eK.\n",
00128 TempNew , phycon.TEMP_LIMIT_LOW );
00129 fprintf(ioQQQ," This calculation is aborting.\n Sorry.\n");
00130 TempNew = phycon.TEMP_LIMIT_LOW*1.0001;
00131 }
00132 else
00133 {
00134
00135 phycon.te = TempNew;
00136 }
00137
00138
00139 tfidle( false );
00140 return;
00141 }
00142
00143 void tfidle(
00144
00145 bool lgForceUpdate)
00146 {
00147 static double tgffused=-1.,
00148 tgffsued2=-1.;
00149 static long int nff_defined=-1;
00150 static long maxion = 0, oldmaxion = 0;
00151 static double ttused = 0.;
00152 static bool lgZLogSet = false;
00153 bool lgGauntF;
00154 long int ion;
00155 long int i,
00156 nelem,
00157 if1,
00158 ipTe,
00159 ret;
00160 double fac,
00161 fanu;
00162
00163 DEBUG_ENTRY( "tfidle()" );
00164
00165
00166 if( lgForceUpdate )
00167 {
00168 ttused = -1.;
00169 tgffused = -1.;
00170 tgffsued2 = -1.;
00171 }
00172
00173
00174 if( dense.eden <= 0. )
00175 {
00176 fprintf( ioQQQ, "tfidle called with a zero or negative electron density,%10.2e\n",
00177 dense.eden );
00178 TotalInsanity();
00179 }
00180
00181
00182 if( phycon.te <= 0. )
00183 {
00184 fprintf( ioQQQ, "tfidle called with a negative electron temperature,%10.2e\n",
00185 phycon.te );
00186 TotalInsanity();
00187 }
00188
00189
00190 if( !lgZLogSet )
00191 {
00192 for( nelem=0; nelem<LIMELM; ++nelem )
00193 {
00194
00195
00196 phycon.sqlogz[nelem] = log10( POW2(nelem+1.) );
00197 }
00198 lgZLogSet = true;
00199 }
00200
00201 if( ! fp_equal( phycon.te, ttused ) )
00202 {
00203 ttused = phycon.te;
00204 thermal.te_update = phycon.te;
00205
00206 phycon.te_eV = phycon.te/EVDEGK;
00207 phycon.te_ryd = phycon.te/TE1RYD;
00208 phycon.te_wn = phycon.te / T1CM;
00209
00210 phycon.tesqrd = POW2(phycon.te);
00211 phycon.sqrte = sqrt(phycon.te);
00212 thermal.halfte = 0.5/phycon.te;
00213 thermal.tsq1 = 1./phycon.tesqrd;
00214 phycon.te32 = phycon.te*phycon.sqrte;
00215 phycon.teinv = 1./phycon.te;
00216
00217 phycon.alogte = log10(phycon.te);
00218 phycon.alnte = log(phycon.te);
00219
00220 phycon.telogn[0] = phycon.alogte;
00221 for( i=1; i < 7; i++ )
00222 {
00223 phycon.telogn[i] = phycon.telogn[i-1]*phycon.telogn[0];
00224 }
00225
00226 phycon.te10 = pow(phycon.te,0.10);
00227 phycon.te20 = phycon.te10 * phycon.te10;
00228 phycon.te30 = phycon.te20 * phycon.te10;
00229 phycon.te40 = phycon.te30 * phycon.te10;
00230 phycon.te70 = phycon.sqrte * phycon.te20;
00231 phycon.te90 = phycon.te70 * phycon.te20;
00232
00233 phycon.te01 = pow(phycon.te,0.01);
00234 phycon.te02 = phycon.te01 * phycon.te01;
00235 phycon.te03 = phycon.te02 * phycon.te01;
00236 phycon.te04 = phycon.te02 * phycon.te02;
00237 phycon.te05 = phycon.te03 * phycon.te02;
00238 phycon.te07 = phycon.te05 * phycon.te02;
00239
00240 phycon.te001 = pow(phycon.te,0.001);
00241 phycon.te002 = phycon.te001 * phycon.te001;
00242 phycon.te003 = phycon.te002 * phycon.te001;
00243 phycon.te004 = phycon.te002 * phycon.te002;
00244 phycon.te005 = phycon.te003 * phycon.te002;
00245 phycon.te007 = phycon.te005 * phycon.te002;
00246
00247 phycon.te0001 = pow(phycon.te ,0.0001);
00248 phycon.te0002 = phycon.te0001 * phycon.te0001;
00249 phycon.te0003 = phycon.te0002 * phycon.te0001;
00250 phycon.te0004 = phycon.te0002 * phycon.te0002;
00251 phycon.te0005 = phycon.te0003 * phycon.te0002;
00252 phycon.te0007 = phycon.te0005 * phycon.te0002;
00253
00254 }
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268 dense.EdenHCorr = dense.eden +
00269
00270 dense.xIonDense[ipHYDROGEN][0]*1.7e-4 * dense.HCorrFac;
00271
00272
00273
00274 dense.EdenHontoHCorr = dense.EdenHCorr;
00275
00276
00277
00278
00279
00280
00281
00282 dense.edensqte = dense.EdenHCorr/phycon.sqrte;
00283 dense.cdsqte = dense.edensqte*COLL_CONST;
00284 dense.SqrtEden = sqrt(dense.eden);
00285
00286
00287
00288
00289 velset();
00290
00291
00292 if( !lgRfieldMalloced )
00293 {
00294 return;
00295 }
00296
00297
00298
00299
00300
00301
00302
00303
00304
00305 if( ! fp_equal(tgffused, phycon.te)
00306 || rfield.ContBoltz[0] <= 0. || nff_defined<rfield.nflux )
00307 {
00308 tgffused = phycon.te;
00309 fac = TE1RYD/phycon.te;
00310 i = 0;
00311 fanu = fac*rfield.anu[i];
00312
00313
00314
00315
00316
00317 while( i < rfield.nupper && fanu < SEXP_LIMIT )
00318 {
00319 rfield.ContBoltz[i] = exp(-fanu);
00320 ++i;
00321
00322 fanu = fac*rfield.anu[i];
00323 }
00324
00325 rfield.ipMaxBolt = i;
00326
00327
00328
00329 for( i=rfield.ipMaxBolt; i < rfield.nupper; i++ )
00330 {
00331 rfield.ContBoltz[i] = 0.;
00332 }
00333 }
00334
00335
00336 tauff();
00337
00338 oldmaxion = maxion;
00339 maxion = 0;
00340
00341
00342 for( nelem = 0; nelem < LIMELM; nelem++ )
00343 {
00344 maxion = MAX2( maxion, dense.IonHigh[nelem] );
00345 }
00346
00347
00348
00349 const double LIM = 0.02;
00350 if( fabs(1.-tgffsued2/phycon.te) > LIM ||
00351
00352
00353
00354 nff_defined<rfield.nflux
00355
00356 || maxion>oldmaxion )
00357 {
00358 static bool lgFirstRunDone = false;
00359 long lowion;
00360
00361 if( fabs(1.-tgffsued2/phycon.te) <= LIM && nff_defined>=rfield.nflux &&
00362 maxion>oldmaxion )
00363 {
00364
00365
00366
00367 lowion = oldmaxion;
00368 }
00369 else
00370 {
00371
00372 lowion = 1;
00373 }
00374
00375
00376 if1 = 1;
00377
00378
00379
00380
00381
00382
00383 nff_defined = rfield.nflux;
00384
00385 ASSERT( if1 >= 0 );
00386
00387 tgffsued2 = phycon.te;
00388 lgGauntF = true;
00389
00390
00391 if( lgGffNotFilled )
00392 {
00393 FillGFF();
00394 }
00395
00396 if( lgFirstRunDone == false )
00397 {
00398 maxion = LIMELM;
00399 lgFirstRunDone = true;
00400 }
00401
00402
00403
00404
00405
00406
00407 ipTe = -1;
00408 for( ion=1; ion<=LIMELM; ion++ )
00409 {
00410
00411 ret = LinterpTable(GauntFF[ion], TeGFF, N_PHOTON_GFF, N_TE_GFF, (realnum)phycon.alogte, GauntFF_T[ion], &ipTe);
00412 if(ret == -1)
00413 {
00414 fprintf(ioQQQ," LinterpTable for GffTable called with te out of bounds \n");
00415 cdEXIT(EXIT_FAILURE);
00416 }
00417 }
00418
00419
00420
00421
00422 ret = LinterpVector(GauntFF_T+lowion, PhoGFF, maxion-lowion+1, N_PHOTON_GFF,
00423 rfield.anulog, rfield.nupper, rfield.gff + lowion);
00424 if(ret == -1)
00425 {
00426 fprintf(ioQQQ," LinterpVector for GffTable called with photon energy out of bounds \n");
00427 cdEXIT(EXIT_FAILURE);
00428 }
00429 }
00430 else
00431 {
00432
00433
00434
00435 if1 = -1;
00436 lgGauntF = false;
00437 }
00438
00439 if( trace.lgTrace && trace.lgTrGant )
00440 {
00441 fprintf( ioQQQ, " tfidle; gaunt factors?" );
00442 fprintf( ioQQQ, "%2c", TorF(lgGauntF) );
00443
00444 fprintf( ioQQQ, "%2f g 1 2=%10.2e%9.2ld flag%2f guv(1,n)%10.2e\n",
00445 rfield.gff[1][0], rfield.gff[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1],
00446 if1, rfield.gff[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]],
00447 rfield.gff[1][rfield.nflux-1] );
00448 }
00449 return;
00450 }
00451
00452
00453 STATIC void tauff(void)
00454 {
00455 realnum fac;
00456
00457
00458 if( !lgOpacMalloced )
00459 return;
00460
00461 DEBUG_ENTRY( "tauff()" );
00462
00463
00464
00465 while( rfield.ipEnergyBremsThin < rfield.nflux &&
00466 opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] >= 1. )
00467 {
00468 ++rfield.ipEnergyBremsThin;
00469 }
00470
00471
00472
00473 if( rfield.ipEnergyBremsThin > 1 && opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] > 0.001 )
00474 {
00475
00476 fac = (1.f - opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1])/(opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] -
00477 opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1]);
00478 fac = MAX2(fac,0.f);
00479 rfield.EnergyBremsThin = rfield.anu[rfield.ipEnergyBremsThin-1] + rfield.widflx[rfield.ipEnergyBremsThin-1]*fac;
00480 }
00481 else
00482 {
00483 rfield.EnergyBremsThin = 0.f;
00484 }
00485
00486
00487 rfield.plsfrq = (realnum)(2.729e-12*sqrt(dense.eden*1.2));
00488
00489 if( rfield.ipPlasma > 0 )
00490 {
00491
00492 while( rfield.plsfrq > rfield.anu[rfield.ipPlasma]+rfield.widflx[rfield.ipPlasma]/2. )
00493 ++rfield.ipPlasma;
00494
00495
00496 while( rfield.ipPlasma>2 && rfield.plsfrq < rfield.anu[rfield.ipPlasma]-rfield.widflx[rfield.ipPlasma]/2. )
00497 --rfield.ipPlasma;
00498 }
00499
00500
00501 rfield.plsfrqmax = MAX2(rfield.plsfrqmax, rfield.plsfrq);
00502
00503
00504 if( rfield.plsfrq > rfield.anu[0] )
00505 {
00506 rfield.lgPlasNu = true;
00507 }
00508
00509
00510
00511 rfield.EnergyBremsThin = MAX2(rfield.plsfrq,rfield.EnergyBremsThin);
00512
00513
00514 while( rfield.ipEnergyBremsThin < rfield.nflux &&
00515 rfield.anu[rfield.ipEnergyBremsThin] <= rfield.EnergyBremsThin )
00516 {
00517 ++rfield.ipEnergyBremsThin;
00518 }
00519 return;
00520 }
00521
00522
00523 STATIC void velset(void)
00524 {
00525 long int nelem;
00526 double turb2;
00527
00528 DEBUG_ENTRY( "velset()" );
00529
00530
00531
00532 turb2 = POW2(DoppVel.TurbVel);
00533
00534
00535
00536
00537 if( DoppVel.DispScale > 0. )
00538 {
00539
00540 turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale );
00541 }
00542
00543
00544
00545 if( wind.windv0 < 0. )
00546 {
00547 turb2 += POW2(wind.windv0);
00548 }
00549
00550
00551
00552 for( nelem=0; nelem < LIMELM; nelem++ )
00553 {
00554 DoppVel.doppler[nelem] =
00555
00556
00557 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/dense.AtomicWeight[nelem]+
00558 turb2);
00559
00560 DoppVel.AveVel[nelem] = sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT*phycon.te/dense.AtomicWeight[nelem]);
00561 }
00562
00563
00564
00565 DoppVel.doppler[LIMELM] =
00566 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
00567 (dense.AtomicWeight[5]+dense.AtomicWeight[7]) + turb2);
00568 DoppVel.AveVel[LIMELM] =
00569 sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
00570 (dense.AtomicWeight[5]+dense.AtomicWeight[7]) + turb2);
00571
00572
00573 DoppVel.doppler[LIMELM+1] =
00574 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
00575 (dense.AtomicWeight[5]*13./12.+dense.AtomicWeight[7]) + turb2);
00576 DoppVel.AveVel[LIMELM+1] =
00577 sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
00578 (dense.AtomicWeight[5]*13./12.+dense.AtomicWeight[7]) + turb2);
00579
00580
00581 DoppVel.doppler[LIMELM+2] =
00582 (realnum)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
00583 (2.*dense.AtomicWeight[0]) + turb2);
00584 DoppVel.AveVel[LIMELM+2] =
00585 sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
00586 (2.*dense.AtomicWeight[0]) + turb2);
00587 return;
00588 }
00589
00590 STATIC void FillGFF( void )
00591 {
00592
00593 long i,i1,i2,i3,j,charge,GffMAGIC = 80321;
00594 double Temp, photon;
00595 bool lgEOL;
00596
00597 DEBUG_ENTRY( "FillGFF()" );
00598
00599 # define chLine_LENGTH 1000
00600 char chLine[chLine_LENGTH];
00601
00602 FILE *ioDATA;
00603
00604 for( i = 0; i < N_TE_GFF; i++ )
00605 {
00606 TeGFF[i] = 0.5f*i;
00607 }
00608
00609 TeGFF[N_TE_GFF-1] += 0.01f;
00610
00611 for( i = 0; i< N_PHOTON_GFF; i++ )
00612 {
00613 PhoGFF[i] = 0.125f*i - 8.f;
00614 }
00615
00616 GauntFF = (realnum***)MALLOC( sizeof(realnum**)*(unsigned)(LIMELM+2) );
00617 for( i = 1; i <= LIMELM; i++ )
00618 {
00619 GauntFF[i] = (realnum**)MALLOC( sizeof(realnum*)*(unsigned)N_PHOTON_GFF );
00620 for( j = 0; j < N_PHOTON_GFF; j++ )
00621 {
00622 GauntFF[i][j] = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_TE_GFF );
00623 }
00624 }
00625
00626 GauntFF_T = (realnum**)MALLOC( sizeof(realnum*)*(unsigned)(LIMELM+2) );
00627 for( i = 1; i <= LIMELM; i++ )
00628 {
00629 GauntFF_T[i] = (realnum*)MALLOC( sizeof(realnum)*(unsigned)N_PHOTON_GFF );
00630 }
00631
00632 if( !rfield.lgCompileGauntFF )
00633 {
00634 if( trace.lgTrace )
00635 fprintf( ioQQQ," FillGFF opening gauntff.dat:");
00636
00637 try
00638 {
00639 ioDATA = open_data( "gauntff.dat", "r" );
00640 }
00641 catch( cloudy_exit )
00642 {
00643 fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
00644 fprintf( ioQQQ, "but the code runs much faster if you compile gauntff.dat!\n");
00645 ioDATA = NULL;
00646 }
00647
00648 if( ioDATA == NULL )
00649 {
00650
00651 for( charge=1; charge<=LIMELM; charge++ )
00652 {
00653 for( i=0; i<N_PHOTON_GFF; i++ )
00654 {
00655 photon = pow((realnum)10.f,PhoGFF[i]);
00656 for(j=0; j<N_TE_GFF; j++)
00657 {
00658 Temp = pow((realnum)10.f,TeGFF[j]);
00659 GauntFF[charge][i][j] = (realnum)cont_gaunt_calc( Temp, (double)charge, photon );
00660 }
00661 }
00662 }
00663 }
00664 else
00665 {
00666
00667 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00668 {
00669 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
00670 cdEXIT(EXIT_FAILURE);
00671 }
00672 i = 1;
00673 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00674 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00675 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00676
00677 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
00678 {
00679 fprintf( ioQQQ,
00680 " FillGFF: the version of gauntff.dat is not the current version.\n" );
00681 fprintf( ioQQQ,
00682 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
00683 GffMAGIC ,
00684 N_PHOTON_GFF,
00685 N_TE_GFF,
00686 i1 , i2 , i3 );
00687 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00688 fprintf( ioQQQ,
00689 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
00690 cdEXIT(EXIT_FAILURE);
00691 }
00692
00693
00694 for( charge = 1; charge <= LIMELM; charge++ )
00695 {
00696 for( i = 0; i<N_PHOTON_GFF; i++ )
00697 {
00698
00699 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00700 {
00701 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
00702 cdEXIT(EXIT_FAILURE);
00703 }
00704
00705 i3 = 1;
00706 i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
00707 i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
00708
00709 if( i1!=charge || i2!=i )
00710 {
00711 fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n");
00712 fprintf( ioQQQ,
00713 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
00714 cdEXIT(EXIT_FAILURE);
00715 }
00716
00717
00718 for(j = 0; j < N_TE_GFF; j++)
00719 {
00720 GauntFF[charge][i][j] = (realnum)FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);
00721
00722 if( lgEOL )
00723 {
00724 fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n");
00725 fprintf( ioQQQ,
00726 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
00727 cdEXIT(EXIT_FAILURE);
00728 }
00729 }
00730 }
00731
00732 }
00733
00734
00735 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00736 {
00737 fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
00738 cdEXIT(EXIT_FAILURE);
00739 }
00740 i = 1;
00741 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00742 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00743 i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00744
00745 if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
00746 {
00747 fprintf( ioQQQ,
00748 " FillGFF: the version of gauntff.dat is not the current version.\n" );
00749 fprintf( ioQQQ,
00750 " FillGFF: I expected to find the numbers %li %li %li and got %li %li %li instead.\n" ,
00751 GffMAGIC ,
00752 N_PHOTON_GFF,
00753 N_TE_GFF,
00754 i1 , i2 , i3 );
00755 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00756 fprintf( ioQQQ,
00757 " FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
00758 cdEXIT(EXIT_FAILURE);
00759 }
00760
00761
00762 fclose( ioDATA );
00763 }
00764 if( trace.lgTrace )
00765 fprintf( ioQQQ," - opened and read ok.\n");
00766 }
00767 else
00768 {
00769
00770
00771 FILE *ioGFF;
00772
00773
00774 ioGFF = open_data( "gauntff.dat" , "w", AS_LOCAL_ONLY );
00775 fprintf(ioGFF,"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
00776 GffMAGIC ,
00777 N_PHOTON_GFF,
00778 N_TE_GFF,
00779 N_PHOTON_GFF,
00780 N_TE_GFF );
00781
00782 for( charge = 1; charge <= LIMELM; charge++ )
00783 {
00784 for( i=0; i < N_PHOTON_GFF; i++ )
00785 {
00786 fprintf(ioGFF, "%li\t%li", charge, i );
00787
00788 for(j = 0; j < N_TE_GFF; j++)
00789 {
00790
00791 Temp = pow((realnum)10.f,TeGFF[j]);
00792 photon = pow((realnum)10.f,PhoGFF[i]);
00793 GauntFF[charge][i][j] = (realnum)cont_gaunt_calc( Temp, (double)charge, photon );
00794 fprintf(ioGFF, "\t%f", GauntFF[charge][i][j] );
00795 }
00796 fprintf(ioGFF, "\n" );
00797 }
00798 }
00799
00800
00801 fprintf(ioGFF,"%li\t%li\t%li\tGFF free-free gaunt factors, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
00802 GffMAGIC ,
00803 N_PHOTON_GFF,
00804 N_TE_GFF,
00805 N_PHOTON_GFF,
00806 N_TE_GFF );
00807
00808 fclose( ioGFF );
00809
00810 fprintf( ioQQQ, "FillGFF: compilation complete, gauntff.dat created.\n" );
00811 }
00812
00813 lgGffNotFilled = false;
00814
00815 {
00816
00817
00818
00819 enum {DEBUG_LOC=false};
00820
00821
00822 if( DEBUG_LOC )
00823 {
00824 double gaunt, error;
00825 double tempsave = phycon.te;
00826 long logu, loggamma2;
00827
00828 for( logu=-4; logu<=4; logu++)
00829 {
00830
00831
00832
00833
00834 for(loggamma2=-4; loggamma2<=4; loggamma2++)
00835 {
00836 double SutherlandGff[9][9]=
00837 { {5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
00838 {4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
00839 {3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
00840 {1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
00841 {0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
00842 {0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
00843 {0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
00844 {0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
00845 {0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};
00846
00847 phycon.te = (TE1RYD/pow(10.,(double)loggamma2));
00848 phycon.alogte = log10(phycon.te);
00849 gaunt = InterpolateGff( 1, pow(10.,(double)(logu-loggamma2)) );
00850 error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
00851
00852 if( error>0.05 )
00853 {
00854 fprintf(ioQQQ," tfidle found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
00855 logu, loggamma2, error );
00856 }
00857 }
00858
00859 }
00860 phycon.te = tempsave;
00861 phycon.alogte = log10(phycon.te);
00862 }
00863 }
00864
00865 return;
00866 }
00867
00868
00869 STATIC realnum InterpolateGff( long charge , double ERyd )
00870 {
00871 double GauntAtLowerPho, GauntAtUpperPho;
00872 static long int ipTe=-1, ipPho=-1;
00873 double gaunt = 0.;
00874 double xmin , xmax;
00875 long i;
00876
00877 DEBUG_ENTRY( "InterpolateGff()" );
00878
00879 if( ipTe < 0 )
00880 {
00881
00882 if( phycon.alogte < TeGFF[0] || phycon.alogte > TeGFF[N_TE_GFF-1] )
00883 {
00884 fprintf(ioQQQ," InterpolateGff called with te out of bounds \n");
00885 cdEXIT(EXIT_FAILURE);
00886 }
00887
00888 for( i=0; i<N_TE_GFF-1; ++i )
00889 {
00890 if( phycon.alogte > TeGFF[i] && phycon.alogte <= TeGFF[i+1] )
00891 {
00892
00893 ipTe = i;
00894 break;
00895 }
00896 }
00897 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) );
00898
00899 }
00900 else if( phycon.alogte < TeGFF[ipTe] )
00901 {
00902
00903 ASSERT( phycon.alogte > TeGFF[0] );
00904
00905 while( phycon.alogte < TeGFF[ipTe] && ipTe > 0)
00906 --ipTe;
00907 }
00908 else if( phycon.alogte > TeGFF[ipTe+1] )
00909 {
00910
00911 ASSERT( phycon.alogte <= TeGFF[N_TE_GFF-1] );
00912
00913 while( phycon.alogte > TeGFF[ipTe+1] && ipTe < N_TE_GFF-1)
00914 ++ipTe;
00915 }
00916
00917 ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1) );
00918
00919
00920 ASSERT( phycon.alogte >= TeGFF[ipTe] && phycon.alogte <= TeGFF[ipTe+1] && ipTe < N_TE_GFF-1 );
00921
00922
00923
00924 if( ipPho < 0 )
00925 {
00926 if( log10(ERyd) < PhoGFF[0] || log10(ERyd) > PhoGFF[N_PHOTON_GFF-1] )
00927 {
00928 fprintf(ioQQQ," InterpolateGff called with photon energy out of bounds \n");
00929 cdEXIT(EXIT_FAILURE);
00930 }
00931 for( i=0; i<N_PHOTON_GFF-1; ++i )
00932 {
00933 if( log10(ERyd) > PhoGFF[i] && log10(ERyd) <= PhoGFF[i+1] )
00934 {
00935 ipPho = i;
00936 break;
00937 }
00938 }
00939 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) );
00940
00941 }
00942 else if( log10(ERyd) < PhoGFF[ipPho] )
00943 {
00944 ASSERT( log10(ERyd) >= PhoGFF[0] );
00945 while( log10(ERyd) < PhoGFF[ipPho] && ipPho > 0)
00946 --ipPho;
00947 }
00948 else if( log10(ERyd) > PhoGFF[ipPho+1] )
00949 {
00950 ASSERT( log10(ERyd) <= PhoGFF[N_PHOTON_GFF-1] );
00951 while( log10(ERyd) > PhoGFF[ipPho+1] && ipPho < N_PHOTON_GFF-1)
00952 ++ipPho;
00953 }
00954 ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1) );
00955 ASSERT( log10(ERyd)>=PhoGFF[ipPho]
00956 && log10(ERyd)<=PhoGFF[ipPho+1] && ipPho<N_PHOTON_GFF-1 );
00957
00958
00959
00960
00961
00962 GauntAtLowerPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) *
00963 (GauntFF[charge][ipPho][ipTe+1] - GauntFF[charge][ipPho][ipTe]) + GauntFF[charge][ipPho][ipTe];
00964
00965 GauntAtUpperPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) *
00966 (GauntFF[charge][ipPho+1][ipTe+1] - GauntFF[charge][ipPho+1][ipTe]) + GauntFF[charge][ipPho+1][ipTe];
00967
00968 gaunt = (log10(ERyd) - PhoGFF[ipPho]) / (PhoGFF[ipPho+1] - PhoGFF[ipPho]) *
00969 (GauntAtUpperPho - GauntAtLowerPho) + GauntAtLowerPho;
00970
00971 xmax = MAX4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1],
00972 GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] );
00973 ASSERT( gaunt <= xmax );
00974
00975 xmin = MIN4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1],
00976 GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] );
00977 ASSERT( gaunt >= xmin );
00978
00979 ASSERT( gaunt > 0. );
00980
00981 return (realnum)gaunt;
00982 }
00983
00984
00985
00986
00987
00988 STATIC int LinterpTable(realnum **t, realnum *v, long int lta, long int ltb, realnum x, realnum *a, long int *pipx)
00989 {
00990 long int ipx=-1;
00991 realnum frac;
00992 long i;
00993
00994 DEBUG_ENTRY( "LinterpTable()" );
00995
00996 if(pipx != NULL)
00997 ipx = *pipx;
00998
00999 fhunt (v,ltb,x,&ipx);
01000 if(pipx != NULL)
01001 *pipx = ipx;
01002
01003 if( ipx == -1 || ipx == ltb )
01004 {
01005 return -1;
01006 }
01007
01008 ASSERT( (ipx >=0) && (ipx < ltb-1) );
01009 ASSERT( x >= v[ipx] && x <= v[ipx+1]);
01010
01011 frac = (x - v[ipx]) / (v[ipx+1] - v[ipx]);
01012 for( i=0; i<lta; i++ )
01013 {
01014 a[i] = frac*t[i][ipx+1]+(1.f-frac)*t[i][ipx];
01015 }
01016
01017 return 0;
01018 }
01019
01020
01021
01022 STATIC int LinterpVector(realnum **t, realnum *v, long lta , long ltb, realnum *yy, long ly, realnum **a)
01023 {
01024 realnum yl, frac;
01025 long i, j, n;
01026
01027 DEBUG_ENTRY( "LinterpVector()" );
01028
01029 if( yy[0] < v[0] || yy[ly-1] > v[ltb-1] )
01030 {
01031 return -1;
01032 }
01033
01034 n = 0;
01035 yl = yy[n];
01036 for(j = 1; j < ltb && n < ly; j++) {
01037 while (yl < v[j] && n < ly) {
01038 frac = (yl-v[j-1])/(v[j]-v[j-1]);
01039 for(i = 0; i < lta; i++)
01040 a[i][n] = frac*t[i][j]+(1.f-frac)*t[i][j-1];
01041 n++;
01042 if(n == ly)
01043 break;
01044 ASSERT( yy[n] >= yy[n-1] );
01045 yl = yy[n];
01046 }
01047 }
01048
01049 return 0;
01050 }
01051 STATIC void fhunt(realnum *xx, long int n, realnum x, long int *j)
01052 {
01053
01054 long int jl, jm, jh, in;
01055 int up;
01056
01057 jl = *j;
01058 up = (xx[n-1] > xx[0]);
01059 if(jl < 0 || jl >= n)
01060 {
01061 jl = -1;
01062 jh = n;
01063 }
01064 else
01065 {
01066 in = 1;
01067 if((x >= xx[jl]) == up)
01068 {
01069 if(jl == n-1)
01070 {
01071 *j = jl;
01072 return;
01073 }
01074 jh = jl + 1;
01075 while ((x >= xx[jh]) == up)
01076 {
01077 jl = jh;
01078 in += in;
01079 jh += in;
01080 if(jh >= n)
01081 {
01082 jh = n;
01083 break;
01084 }
01085 }
01086 }
01087 else
01088 {
01089 if(jl == 0)
01090 {
01091 *j = -1;
01092 return;
01093 }
01094 jh = jl--;
01095 while ((x < xx[jl]) == up)
01096 {
01097 jh = jl;
01098 in += in;
01099 jl -= in;
01100 if(jl <= 0)
01101 {
01102 jl = 0;
01103 break;
01104 }
01105 }
01106 }
01107 }
01108 while (jh-jl != 1)
01109 {
01110 jm = (jh+jl)/2;
01111 if((x > xx[jm]) == up)
01112 jl = jm;
01113 else
01114 jh = jm;
01115 }
01116 *j = jl;
01117 return;
01118 }
01119