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