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