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