00001
00002
00003
00004
00005
00006
00007 #include "cddefines.h"
00008 #include "cddrive.h"
00009 #include "physconst.h"
00010 #include "optimize.h"
00011 #include "rfield.h"
00012 #include "trace.h"
00013 #include "radius.h"
00014 #include "input.h"
00015 #include "stars.h"
00016 #include "lines.h"
00017 #include "prt.h"
00018 #include "parser.h"
00019 #include "save.h"
00020 #include "thirdparty.h"
00021 #include "continuum.h"
00022
00023 #if defined(__HP_aCC)
00024 # pragma OPT_LEVEL 1
00025 #endif
00026
00027
00028 STATIC void ReadTable(const char * fnam);
00029
00030 static string chLINE_LIST;
00031
00032
00033
00034 static const int NCRAB = 10;
00035 static double tnucrb[NCRAB],
00036 fnucrb[NCRAB];
00037
00038
00039 static const int NRUBIN = 56;
00040 static double tnurbn[NRUBIN] = {1.05E-08,1.05E-07,1.05E-06,1.04E-05,1.00E-04,1.00E-03,1.00E-02,3.01E-02,1.00E-01,
00041 1.50E-01,2.50E-01,4.01E-01,6.01E-01,9.8E-01,9.96E-01,1.00E+00,1.02445,1.07266,1.12563,1.18411,1.23881,
00042 1.29328,1.35881,1.42463,1.48981,1.55326,1.6166,1.68845,1.76698,1.8019,1.808,1.84567,1.9317,2.04891,2.14533,
00043 2.19702,2.27941,2.37438,2.43137,2.49798,2.56113,2.59762,2.66597,2.80543,2.95069,3.02911,3.11182,3.22178,
00044 3.3155,3.42768,3.50678,3.56157,3.61811,3.75211,3.89643,4.},
00045 fnurbn[NRUBIN] = {1.56E-20,1.55E-17,1.54E-14,1.53E-11,1.35E-08,1.34E-05,1.35E-02,3.62E-01,1.27E+01,
00046 3.90E+01,1.48E+02,4.52E+02,1.02E+03,2.27E+03,8.69E+02,8.04E+02,6.58E+02,6.13E+02,6.49E+02,6.06E+02,
00047 6.30E+02,5.53E+02,5.55E+02,5.24E+02,4.86E+02,4.49E+02,4.42E+02,3.82E+02,3.91E+02,2.91E+02,2.61E+02,
00048 1.32E+02,1.20E+02,1.16E+02,9.56E+01,9.94E+01,9.10E+01,4.85E+01,7.53E+01,9.53E+01,8.52E+01,5.76E+01,
00049 6.72E+01,5.20E+01,8.09E+00,3.75E+00,5.57E+00,3.80E+00,2.73E+00,2.22E+00,3.16E-01,1.26E-01,1.39E-01,
00050 6.15E-02,3.22E-02,7.98E-03};
00051
00052
00053 static const int NCFL = 40;
00054 static double cfle[NCFL],
00055 cflf[NCFL];
00056
00057
00058 static const int NKN120 = 11;
00059 static double tnuakn[NKN120],
00060 fnuakn[NKN120];
00061
00062
00063 static const int NISM = 23;
00064 static double tnuism[NISM],
00065 fnuism[NISM];
00066
00067
00068
00069
00070 static const int NHM96 = 14;
00071
00072 static const double tnuHM96[NHM96]={-8,-1.722735683,-0.351545683,-0.222905683,-0.133385683,
00073
00074 -0.127655683,-0.004575683,0.297544317,0.476753,0.476756,0.588704317,
00075 0.661374317,1.500814317,2.245164317};
00076
00077
00078 static const double fnuHM96[NHM96]={-32.53342863,-19.9789,-20.4204,-20.4443,-20.5756,-20.7546,
00079 -21.2796,-21.6256,-21.8404,-21.4823,-22.2102,-22.9263,-23.32,-24.2865};
00080
00081
00082 static const int NAGN = 8;
00083 static double tnuagn[NAGN],
00084 tslagn[NAGN];
00085
00086
00087 static const int NDRAINE = 15;
00088 static double tnudrn[NDRAINE] , tsldrn[NDRAINE];
00089
00090
00091 STATIC void ZeroContin(void);
00092
00093
00094
00095
00096
00097
00098 STATIC void resetBltin( double *tnu , double *fluxlog , bool lgLog )
00099 {
00100
00101
00102 const double RESETFACTOR = 0.98;
00103 double power;
00104
00105 ASSERT( rfield.emm > 0. );
00106
00107 if( lgLog )
00108 {
00109
00110
00111 power = (fluxlog[1] - fluxlog[0] ) / log10( tnu[1]/tnu[0] );
00112
00113 tnu[0] = rfield.emm*RESETFACTOR;
00114 fluxlog[0] = fluxlog[1] + power * log10( tnu[0] /tnu[1] );
00115 }
00116 else
00117 {
00118
00119
00120 power = log10( fluxlog[1]/fluxlog[0]) / log10( tnu[1]/tnu[0] );
00121
00122 tnu[0] = rfield.emm*RESETFACTOR;
00123 fluxlog[0] = log10(fluxlog[1]) + power * log10( tnu[0] /tnu[1] );
00124
00125 fluxlog[0] = pow(10. , fluxlog[0]);
00126 }
00127
00128 return;
00129 }
00130
00131
00132
00133 STATIC void read_hm05( const char chFile[] , double **tnuHM , double **fnuHM ,
00134 long int *nhm , double redshift )
00135 {
00136
00137 FILE *ioFILE;
00138 double *table_redshifts = NULL,
00139 *table_wl = NULL ,
00140 **table_fn=NULL,
00141 frac_hi;
00142 char chLine[1000];
00143 long int nRedshift , i , n , nz , ipLo , ipHi;
00144 bool lgEOL;
00145
00146 DEBUG_ENTRY( "read_hm05()" );
00147
00148 ioFILE = open_data( chFile, "r", AS_LOCAL_DATA );
00149
00150
00151 *nhm = 0;
00152
00153 nRedshift = -1;
00154
00155 while( read_whole_line( chLine , (int)sizeof(chLine) , ioFILE ) != NULL )
00156 {
00157
00158
00159 if( chLine[0] != '#')
00160 {
00161 ++*nhm;
00162
00163 if( *nhm==1 )
00164 {
00165 long mag_read;
00166
00167 i = 1;
00168 mag_read = (long)FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
00169 if( mag_read != 50808 )
00170 {
00171 fprintf(ioQQQ,
00172 " Magic number in Haardt & Madau file is not correct, I expected 50808 and found %li\n",
00173 mag_read );
00174 cdEXIT(EXIT_FAILURE);
00175 }
00176 }
00177
00178 else if( *nhm==2 )
00179 {
00180 double a;
00181 i = 2;
00182 lgEOL = false;
00183 nRedshift = 0;
00184 while( !lgEOL )
00185 {
00186 ++nRedshift;
00187 a = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
00188 ASSERT( a >= 0. );
00189 }
00190
00191 --nRedshift;
00192
00193 --nRedshift;
00194
00195
00196 table_redshifts = (double*)MALLOC( (size_t)nRedshift*sizeof(double) );
00197
00198 i = 2;
00199 for( n=0; n<nRedshift; ++n )
00200 {
00201 table_redshifts[n] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
00202
00203 }
00204 if( lgEOL )
00205 TotalInsanity();
00206 }
00207 }
00208 }
00209
00210
00211 table_wl = (double*)MALLOC( (size_t)*nhm*sizeof(double) );
00212
00213 table_fn = (double**)MALLOC( (size_t)nRedshift*sizeof(double*) );
00214 for(n=0; n<nRedshift; ++n )
00215 {
00216 table_fn[n] = (double*)MALLOC( (size_t)(*nhm)*sizeof(double) );
00217 }
00218
00219
00220 if( fseek( ioFILE , 0 , SEEK_SET ) != 0 )
00221 {
00222 fprintf( ioQQQ, " read_hm05 could not rewind HM05 date file.\n");
00223 cdEXIT(EXIT_FAILURE);
00224 }
00225 n = 0;
00226 *nhm = 0;
00227 while( read_whole_line( chLine , (int)sizeof(chLine) , ioFILE ) != NULL )
00228 {
00229
00230
00231 if( chLine[0] != '#')
00232 {
00233
00234
00235 ++n;
00236 if( n>2 )
00237 {
00238 i = 1;
00239 table_wl[*nhm] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
00240 if( lgEOL )
00241 TotalInsanity();
00242 for( nz=0; nz<nRedshift; ++nz )
00243 {
00244
00245 table_fn[nz][*nhm] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL);
00246 }
00247 ++*nhm;
00248 }
00249 }
00250 }
00251
00252 fclose( ioFILE );
00253
00254 {
00255
00256 enum {DEBUG_LOC=false};
00257 if( DEBUG_LOC)
00258 {
00259 fprintf(ioQQQ,"wavelength/z");
00260 for(nz=0; nz<nRedshift; ++nz )
00261 fprintf(ioQQQ,"\t%.3f", table_redshifts[nz] );
00262 fprintf(ioQQQ,"\n");
00263 for( i=0; i<*nhm; ++i )
00264 {
00265 fprintf(ioQQQ,"%.3e", table_wl[i] );
00266 for( nz=0; nz<nRedshift; ++nz )
00267 fprintf(ioQQQ,"\t%.3e", table_fn[nz][i] );
00268 fprintf(ioQQQ,"\n");
00269 }
00270 }
00271 }
00272
00273
00274 assert( table_redshifts!=NULL );
00275
00276
00277 if( redshift < table_redshifts[0] ||
00278 redshift > table_redshifts[nRedshift-1] )
00279 {
00280 fprintf(ioQQQ," The redshift requested on table HM05 is %g but is not within the range of the table, which goes from %g to %g.\n",
00281 redshift, table_redshifts[0] , table_redshifts[nRedshift-1] );
00282 fprintf(ioQQQ," Sorry.\n");
00283 cdEXIT(EXIT_FAILURE);
00284 }
00285
00286 ipLo = -1;
00287 ipHi = -1;
00288
00289 for( nz=0; nz<nRedshift-1; ++nz )
00290 {
00291 if( redshift >= table_redshifts[nz] &&
00292 redshift <= table_redshifts[nz+1] )
00293 {
00294 ipLo = nz;
00295 ipHi = nz+1;
00296 break;
00297 }
00298 }
00299 ASSERT( ipLo>=0 && ipHi >=0 );
00300
00301
00302
00303 for( i=0; i<*nhm-1; ++i )
00304 {
00305 if( table_wl[i]>=table_wl[i+1] )
00306 {
00307 fprintf(ioQQQ," The wavelengths in the table HM05 data table are not in increasing order. This is required.\n");
00308 fprintf(ioQQQ," Sorry.\n");
00309 cdEXIT(EXIT_FAILURE);
00310 }
00311 }
00312
00313
00314 *fnuHM = (double *)MALLOC( (size_t)(*nhm)*sizeof(double ) );
00315 *tnuHM = (double *)MALLOC( (size_t)(*nhm)*sizeof(double ) );
00316
00317
00318
00319
00320 frac_hi = (redshift-table_redshifts[ipLo]) / (table_redshifts[ipHi]-table_redshifts[ipLo]);
00321 for( i=0; i<*nhm; ++i )
00322 {
00323
00324
00325
00326 (*tnuHM)[*nhm-1-i] = RYDLAM / table_wl[i];
00327
00328 (*fnuHM)[*nhm-1-i] = table_fn[ipLo][i]*(1.-frac_hi) + table_fn[ipHi][i]*frac_hi;
00329
00330 }
00331
00332 for( n=0; n<nRedshift; ++n )
00333 free( table_fn[n] );
00334 free( table_fn );
00335 free( table_wl );
00336 free( table_redshifts );
00337 return;
00338 }
00339
00340 void ParseTable(Parser &p)
00341 {
00342 char chFile[FILENAME_PATH_LENGTH_2];
00343
00344 IntMode imode = IM_ILLEGAL_MODE;
00345 bool lgHit,
00346 lgLogSet;
00347 static bool lgCalled=false;
00348
00349 long int i,
00350 j,
00351 nstar,
00352 ipNORM;
00353
00354 double alpha,
00355 brakmm,
00356 brakxr,
00357 ConBreak,
00358 fac,
00359 scale,
00360 slopir,
00361 slopxr;
00362
00363 bool lgNoContinuum = false,
00364 lgQuoteFound;
00365
00366 DEBUG_ENTRY( "ParseTable()" );
00367
00368
00369 if( !lgCalled )
00370 {
00371 ZeroContin();
00372 lgCalled = true;
00373 }
00374
00375 if( rfield.nShape >= LIMSPC )
00376 {
00377 fprintf( ioQQQ, " %ld is too many spectra entered. Increase LIMSPC\n Sorry.\n",
00378 rfield.nShape );
00379 cdEXIT(EXIT_FAILURE);
00380 }
00381
00382
00383
00384
00385 lgQuoteFound = true;
00386 if( p.GetQuote( chFile , false ) )
00387 lgQuoteFound = false;
00388
00389
00390 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
00391
00392 if( !rfield.lgContMalloc[rfield.nShape] )
00393 {
00394 rfield.tNu[rfield.nShape].resize( NCELL );
00395 rfield.tslop[rfield.nShape].resize( NCELL );
00396 rfield.tFluxLog[rfield.nShape].resize( NCELL );
00397 rfield.ncont[rfield.nShape] = 0;
00398 rfield.lgContMalloc[rfield.nShape] = true;
00399 }
00400
00401
00402 if( p.nMatch(" AGN") )
00403 {
00404
00405 ASSERT( NAGN < NCELL);
00406 for( i=0; i < NAGN; i++ )
00407 {
00408 rfield.tNu[rfield.nShape][i].set(tnuagn[i]);
00409 rfield.tslop[rfield.nShape][i] = (realnum)tslagn[i];
00410 }
00411 rfield.ncont[rfield.nShape] = NAGN;
00412
00413
00414 if( p.nMatch("BREA") )
00415 {
00416 ConBreak = p.FFmtRead();
00417
00418 if( p.lgEOL() )
00419 {
00420
00421 if( p.nMatch(" NO ") )
00422 {
00423 ConBreak = rfield.emm*1.01f;
00424 }
00425 else
00426 {
00427 fprintf( ioQQQ, " There must be a number for the break.\n Sorry.\n" );
00428 cdEXIT(EXIT_FAILURE);
00429 }
00430 }
00431
00432 if( ConBreak == 0. )
00433 {
00434 fprintf( ioQQQ, " The break must be greater than 0.2 Ryd.\n Sorry.\n" );
00435 cdEXIT(EXIT_FAILURE);
00436 }
00437
00438 if( p.nMatch("MICR") )
00439 {
00440
00441 ConBreak = 0.0912/ConBreak;
00442 }
00443
00444 if( ConBreak < 0. )
00445 {
00446
00447 ConBreak = pow(10.,ConBreak);
00448 }
00449
00450 else if( ConBreak == 0. )
00451 {
00452 fprintf( ioQQQ, " An energy of 0 is not allowed.\n Sorry.\n" );
00453 cdEXIT(EXIT_FAILURE);
00454 }
00455
00456 if( ConBreak >= rfield.tNu[rfield.nShape][2].Ryd() )
00457 {
00458 fprintf( ioQQQ, " The energy of the break cannot be greater than%10.2e Ryd.\n Sorry.\n",
00459 rfield.tNu[rfield.nShape][2].Ryd() );
00460 cdEXIT(EXIT_FAILURE);
00461 }
00462
00463 else if( ConBreak <= rfield.tNu[rfield.nShape][0].Ryd() )
00464 {
00465 fprintf( ioQQQ, " The energy of the break cannot be less than%10.2e Ryd.\n Sorry.\n",
00466 rfield.tNu[rfield.nShape][0].Ryd() );
00467 cdEXIT(EXIT_FAILURE);
00468 }
00469
00470 rfield.tNu[rfield.nShape][1].set(ConBreak);
00471
00472 rfield.tslop[rfield.nShape][1] =
00473 (realnum)(rfield.tslop[rfield.nShape][2] +
00474 log10(rfield.tNu[rfield.nShape][2].Ryd()/rfield.tNu[rfield.nShape][1].Ryd()));
00475
00476 rfield.tslop[rfield.nShape][0] =
00477 (realnum)(rfield.tslop[rfield.nShape][1] -
00478 2.5*log10(rfield.tNu[rfield.nShape][1].Ryd()/rfield.tNu[rfield.nShape][0].Ryd()));
00479 }
00480 }
00481
00482 else if( p.nMatchErase("AKN1") )
00483 {
00484
00485 ASSERT( NKN120 < NCELL );
00486 for( i=0; i < NKN120; i++ )
00487 {
00488 rfield.tNu[rfield.nShape][i].set(tnuakn[i]);
00489 rfield.tslop[rfield.nShape][i] = (realnum)log10(fnuakn[i]);
00490 }
00491 rfield.ncont[rfield.nShape] = NKN120;
00492 }
00493
00494 else if( p.nMatch("CRAB") )
00495 {
00496 if( p.nMatch("DAVIDSON") )
00497 {
00498 ASSERT( NCRAB < NCELL );
00499
00500 for( i=0; i < NCRAB; i++ )
00501 {
00502 rfield.tNu[rfield.nShape][i].set(tnucrb[i]);
00503 rfield.tslop[rfield.nShape][i] = (realnum)log10(fnucrb[i]);
00504 }
00505 rfield.ncont[rfield.nShape] = NCRAB;
00506 }
00507 else
00508 {
00509
00510
00511
00512 const int NCRAB08 = 14;
00513 ASSERT( NCRAB08 < NCELL );
00514 static const double tnucrbHz08[NCRAB08] = {
00515 10.0114,
00516 12.0499,
00517 13.8657,
00518 14.8718,
00519 15.6126,
00520 16.3815,
00521 18.0873,
00522 18.6467,
00523 19.6535,
00524 20.9825,
00525 21.8082,
00526 22.5783,
00527 23.293,
00528 23.644 };
00529 static const double crbnuLnu08[NCRAB08] = {
00530 34.4381,
00531 35.855,
00532 36.7703,
00533 37.0142,
00534 37.0441,
00535 37.0297,
00536 36.8609,
00537 36.7579,
00538 36.5962,
00539 36.0585,
00540 35.5279,
00541 34.8277,
00542 33.862,
00543 33.0141 };
00544 static double tnu[NCRAB08] , tflux[NCRAB08];
00545 for( i=0; i < NCRAB08; i++ )
00546 {
00547
00548 tnu[i] = pow(10. , tnucrbHz08[i] ) / FR1RYD;
00549
00550 tflux[i] = crbnuLnu08[i] - tnucrbHz08[i];
00551 }
00552 resetBltin( tnu , tflux , false );
00553 for( i=0; i < NCRAB08; i++ )
00554 {
00555 rfield.tNu[rfield.nShape][i].set(tnu[i]);
00556
00557 rfield.tslop[rfield.nShape][i] = (realnum)tflux[i];
00558 }
00559 rfield.ncont[rfield.nShape] = NCRAB08;
00560 }
00561 }
00562
00563 else if( p.nMatch("COOL") )
00564 {
00565 ASSERT( NCFL < NCELL );
00566
00567 for( i=0; i < NCFL; i++ )
00568 {
00569 rfield.tNu[rfield.nShape][i].set(cfle[i]);
00570 rfield.tslop[rfield.nShape][i] = (realnum)cflf[i];
00571 }
00572 rfield.ncont[rfield.nShape] = NCFL;
00573 }
00574
00575 else if( p.nMatchErase("HM96") )
00576 {
00577
00578
00579
00580
00581 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
00582
00583 ASSERT( NHM96 < NCELL );
00584
00585
00586 for( j=0; j < NHM96; j++ )
00587 {
00588
00589 rfield.tNu[rfield.nShape][j].set(pow( 10. , tnuHM96[j] ));
00590 rfield.tslop[rfield.nShape][j] = (realnum)fnuHM96[j];
00591 }
00592 rfield.ncont[rfield.nShape] = NHM96;
00593
00594
00595
00596 scale = p.FFmtRead();
00597 if( scale > 0. )
00598 scale = log10(scale);
00599
00600
00601 if( p.m_nqh >= LIMSPC )
00602 {
00603 fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
00604 p.m_nqh );
00605 cdEXIT(EXIT_FAILURE);
00606 }
00607
00608
00609
00610
00611 if( rfield.nShape != p.m_nqh )
00612 {
00613 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
00614 fprintf( ioQQQ, " Sorry.\n" );
00615 cdEXIT(EXIT_FAILURE);
00616 }
00617
00618 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00619 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
00620
00621
00622 rfield.lgBeamed[p.m_nqh] = false;
00623 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
00624
00625
00626
00627
00628 rfield.range[p.m_nqh][0] = pow(10. , tnuHM96[2] )*1.0001;
00629
00630
00631 rfield.totpow[p.m_nqh] = (fnuHM96[2] + log10(PI4) + scale);
00632
00633
00634
00635 if( !radius.lgRadiusKnown )
00636 {
00637 radius.Radius = pow(10.,radius.rdfalt);
00638 }
00639 ++p.m_nqh;
00640
00641
00642
00643 if( !radius.lgRadiusKnown )
00644 {
00645 radius.Radius = pow(10.,radius.rdfalt);
00646 }
00647 }
00648
00649 else if( p.nMatchErase("HM05") )
00650 {
00651 double *tnuHM , *fnuHM;
00652 double redshift;
00653 long int nhm;
00654
00655
00656
00657
00658
00659 strcpy( rfield.chSpType[rfield.nShape], "INTER" );
00660 if( p.nMatch("QUAS") )
00661 {
00662
00663 strcpy( chFile , "haardt_madau_quasar.dat" );
00664 }
00665 else
00666 {
00667
00668 strcpy( chFile , "haardt_madau_galaxy.dat" );
00669 }
00670
00671
00672
00673 redshift = p.FFmtRead();
00674 if( p.lgEOL() )
00675 {
00676 fprintf(ioQQQ," The redshift MUST be specified on the table HM05 command. I did not find one.\n Sorry.\n");
00677 cdEXIT(EXIT_FAILURE);
00678 }
00679
00680
00681
00682 read_hm05( chFile , &tnuHM , &fnuHM , &nhm , redshift );
00683
00684 ipNORM = -1;
00685 for( j=0; j < nhm-1; j++ )
00686 {
00687
00688
00689 if( tnuHM[j]<=1. && 1. <= tnuHM[j+1] )
00690 {
00691 ipNORM = j;
00692 break;
00693 }
00694 }
00695 ASSERT( ipNORM>=0 );
00696
00697
00698
00699
00700 scale = p.FFmtRead();
00701 if( scale > 0. )
00702 scale = log10(scale);
00703
00704
00705
00706 rfield.range[p.m_nqh][0] = tnuHM[ipNORM];
00707
00708
00709 rfield.totpow[p.m_nqh] = log10(fnuHM[ipNORM]) + log10(PI4) + scale;
00710
00711
00712 if( p.m_nqh >= LIMSPC )
00713 {
00714 fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
00715 p.m_nqh );
00716 cdEXIT(EXIT_FAILURE);
00717 }
00718
00719
00720
00721
00722 if( rfield.nShape != p.m_nqh )
00723 {
00724 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
00725 fprintf( ioQQQ, " Sorry.\n" );
00726 cdEXIT(EXIT_FAILURE);
00727 }
00728
00729 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00730 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
00731
00732
00733 rfield.lgBeamed[p.m_nqh] = false;
00734 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
00735
00736
00737
00738 scale = SDIV( fnuHM[ipNORM] );
00739
00740
00741
00742
00743 rfield.tNu[rfield.nShape][0].set(rfield.emm);
00744 rfield.tslop[rfield.nShape][0] =
00745 (realnum)log10(fnuHM[0]*pow((double)(rfield.emm/tnuHM[0]),2.)/scale);
00746
00747
00748
00749 for( j=0; j < nhm; j++ )
00750 {
00751
00752 rfield.tNu[rfield.nShape][j+1].set(tnuHM[j]);
00753 rfield.tslop[rfield.nShape][j+1] = (realnum)log10(fnuHM[j]/scale);
00754
00755
00756
00757 }
00758 ++nhm;
00759 rfield.ncont[rfield.nShape] = nhm;
00760
00761
00762 for( j=0; j < nhm-1; j++ )
00763 ASSERT( rfield.tNu[rfield.nShape][j].Ryd() < rfield.tNu[rfield.nShape][j+1].Ryd() );
00764
00765
00766
00767 if( !radius.lgRadiusKnown )
00768 {
00769 radius.Radius = pow(10.,radius.rdfalt);
00770 }
00771 ++p.m_nqh;
00772
00773
00774
00775 if( !radius.lgRadiusKnown )
00776 {
00777 radius.Radius = pow(10.,radius.rdfalt);
00778 }
00779 free( tnuHM );
00780 free( fnuHM );
00781
00782 }
00783 else if( p.nMatch(" ISM") )
00784 {
00785 ASSERT( NISM < NCELL );
00786
00787
00788
00789 rfield.tNu[rfield.nShape][0].set(6.);
00790 rfield.tslop[rfield.nShape][0] = -21.21f - 6.f;
00791 for( i=6; i < NISM; i++ )
00792 {
00793 rfield.tNu[rfield.nShape][i-5].set(tnuism[i]);
00794 rfield.tslop[rfield.nShape][i-5] = (realnum)(fnuism[i] -
00795 tnuism[i]);
00796 }
00797
00798 rfield.ncont[rfield.nShape] = NISM -5;
00799
00800
00801
00802
00803
00804 scale = p.FFmtRead();
00805 if( scale > 0. && !p.nMatch(" LOG") )
00806 scale = log10(scale);
00807
00808
00809 if( p.m_nqh >= LIMSPC )
00810 {
00811 fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
00812 p.m_nqh );
00813 cdEXIT(EXIT_FAILURE);
00814 }
00815
00816
00817
00818
00819 if( rfield.nShape != p.m_nqh )
00820 {
00821 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
00822 fprintf( ioQQQ, " Sorry.\n" );
00823 cdEXIT(EXIT_FAILURE);
00824 }
00825
00826 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00827 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
00828
00829
00830 rfield.lgBeamed[p.m_nqh] = false;
00831 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
00832
00833
00834
00835
00836 rfield.range[p.m_nqh][0] = HIONPOT;
00837
00838
00839 rfield.totpow[p.m_nqh] = (-18.517 + scale);
00840
00841
00842
00843 if( !radius.lgRadiusKnown )
00844 {
00845 radius.Radius = pow(10.,radius.rdfalt);
00846 }
00847 ++p.m_nqh;
00848
00849 if( optimize.lgVarOn )
00850 {
00851 optimize.nvarxt[optimize.nparm] = 1;
00852 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE ISM LOG %f");
00853
00854 optimize.nvfpnt[optimize.nparm] = input.nRead;
00855
00856 optimize.vparm[0][optimize.nparm] = (realnum)scale;
00857 optimize.vincr[optimize.nparm] = 0.2f;
00858 ++optimize.nparm;
00859 }
00860 }
00861 else if( p.nMatch("DRAI") )
00862 {
00863 ASSERT( NDRAINE < NCELL );
00864 rfield.lgMustBlockHIon = true;
00865
00866
00867 for( i=0; i < NDRAINE; i++ )
00868 {
00869 rfield.tNu[rfield.nShape][i].set(tnudrn[i]);
00870 rfield.tslop[rfield.nShape][i] = (realnum)tsldrn[i];
00871 }
00872
00873 rfield.ncont[rfield.nShape] = NDRAINE;
00874
00875
00876
00877
00878 scale = p.FFmtRead();
00879 if( scale > 0. && !p.nMatch(" LOG") )
00880 scale = log10(scale);
00881
00882
00883 if( p.m_nqh >= LIMSPC )
00884 {
00885 fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n",
00886 p.m_nqh );
00887 cdEXIT(EXIT_FAILURE);
00888 }
00889
00890
00891
00892
00893 if( rfield.nShape != p.m_nqh )
00894 {
00895 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
00896 fprintf( ioQQQ, " Sorry.\n" );
00897 cdEXIT(EXIT_FAILURE);
00898 }
00899
00900 strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
00901 strcpy( rfield.chSpNorm[p.m_nqh], "FLUX" );
00902
00903
00904 rfield.lgBeamed[p.m_nqh] = false;
00905 rfield.Illumination[p.m_nqh] = Illuminate::SYMMETRIC;
00906
00907
00908
00909
00910 rfield.range[p.m_nqh][0] = tnudrn[0]*1.01;
00911
00912
00913 rfield.totpow[p.m_nqh] = tsldrn[0] + scale;
00914
00915 if( optimize.lgVarOn )
00916 {
00917 optimize.nvarxt[optimize.nparm] = 1;
00918 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE DRAINE LOG %f");
00919
00920 optimize.nvfpnt[optimize.nparm] = input.nRead;
00921
00922 optimize.vparm[0][optimize.nparm] = (realnum)scale;
00923 optimize.vincr[optimize.nparm] = 0.2f;
00924 ++optimize.nparm;
00925 }
00926
00927
00928
00929 if( !radius.lgRadiusKnown )
00930 {
00931 radius.Radius = pow(10.,radius.rdfalt);
00932 }
00933 ++p.m_nqh;
00934 }
00935
00936 else if( p.nMatch("LINE") )
00937 {
00938
00939
00940
00941
00942
00943 lgNoContinuum = true;
00944
00945 if( chLINE_LIST.size() > 0 )
00946 {
00947 fprintf(ioQQQ," sorry, only one table line per input stream\n");
00948 cdEXIT(EXIT_FAILURE);
00949 }
00950
00951
00952
00953 if( lgQuoteFound && strlen(chFile) > 0 )
00954 chLINE_LIST = chFile;
00955 else
00956 chLINE_LIST = "LineList_BLR.dat";
00957
00958
00959
00960 rfield.ncont[rfield.nShape] = -9;
00961
00962
00963 FILE* ioData = open_data( chLINE_LIST.c_str(), "r", AS_LOCAL_DATA_TRY );
00964 if( ioData == NULL )
00965 {
00966
00967 fprintf(ioQQQ,"\n DISASTER PROBLEM ParseTable could not find "
00968 "line list file %s\n", chLINE_LIST.c_str() );
00969 fprintf(ioQQQ," Please check the spelling of the file name and that it "
00970 "is in either the local or data directory.\n\n");
00971 cdEXIT(EXIT_FAILURE);
00972 }
00973 else
00974 {
00975 fclose(ioData);
00976 }
00977
00978 }
00979
00980 else if( p.nMatch("POWE") )
00981 {
00982
00983
00984 alpha = p.FFmtRead();
00985
00986
00987 if( p.lgEOL() )
00988 alpha = -1.;
00989
00990
00991 rfield.tNu[rfield.nShape][0].set(rfield.emm);
00992
00993 rfield.tslop[rfield.nShape][0] = -5.f;
00994
00995 lgLogSet = false;
00996
00997
00998 brakmm = p.FFmtRead();
00999
01000
01001 if( p.lgEOL() )
01002 {
01003 lgLogSet = true;
01004 brakmm = 9.115e-3;
01005 }
01006
01007 else if( brakmm == 0. )
01008 {
01009
01010
01011
01012 lgLogSet = false;
01013 brakmm = rfield.tNu[rfield.nShape][0].Ryd()*1.001;
01014 }
01015
01016 else if( brakmm < 0. )
01017 {
01018
01019 lgLogSet = true;
01020 brakmm = pow(10.,brakmm);
01021 }
01022
01023
01024 if( p.nMatch("MICR") )
01025 brakmm = RYDLAM / (1e4*brakmm);
01026
01027 rfield.tNu[rfield.nShape][1].set(brakmm);
01028
01029
01030 if( brakmm > 1. )
01031 fprintf(ioQQQ,
01032 " Check the order of parameters on this table power law command - the low-energy break of %f Ryd seems high to me.\n",
01033 brakmm );
01034
01035
01036
01037
01038
01039
01040 slopir = 5./2.;
01041
01042
01043
01044 rfield.tslop[rfield.nShape][1] =
01045 (realnum)(rfield.tslop[rfield.nShape][0] +
01046 slopir*log10(rfield.tNu[rfield.nShape][1].Ryd()/rfield.tNu[rfield.nShape][0].Ryd()));
01047
01048
01049 rfield.tNu[rfield.nShape][3].set(rfield.egamry);
01050
01051
01052 brakxr = p.FFmtRead();
01053
01054
01055 if( p.lgEOL() )
01056 {
01057 brakxr = 3676.;
01058 }
01059
01060 else if( brakxr == 0. )
01061 {
01062 brakxr = rfield.tNu[rfield.nShape][3].Ryd()/1.001;
01063 }
01064
01065 else if( lgLogSet )
01066 {
01067
01068 brakxr = pow(10.,brakxr);
01069 }
01070
01071
01072 rfield.tNu[rfield.nShape][2].set(brakxr);
01073
01074
01075
01076 if( brakmm >= brakxr )
01077 {
01078 fprintf( ioQQQ, " HELP!! The lower energy for the power law, %f, is greater than the upper energy, %f. This is not possible.\n Sorry.\n",
01079 brakmm , brakxr );
01080 cdEXIT(EXIT_FAILURE);
01081 }
01082
01083
01084 rfield.tslop[rfield.nShape][2] =
01085 (realnum)(rfield.tslop[rfield.nShape][1] +
01086 alpha*log10(rfield.tNu[rfield.nShape][2].Ryd()/rfield.tNu[rfield.nShape][1].Ryd()));
01087
01088
01089 slopxr = -2.;
01090
01091 rfield.tslop[rfield.nShape][3] =
01092 (realnum)(rfield.tslop[rfield.nShape][2] +
01093 slopxr*log10(rfield.tNu[rfield.nShape][3].Ryd()/rfield.tNu[rfield.nShape][2].Ryd()));
01094
01095
01096 rfield.ncont[rfield.nShape] = 4;
01097 }
01098
01099 else if( p.nMatch("READ") )
01100 {
01101
01102
01103
01104 if( !lgQuoteFound )
01105 {
01106 fprintf( ioQQQ, " Name of file must appear on TABLE READ.\n");
01107 cdEXIT(EXIT_FAILURE);
01108 }
01109
01110
01111 strcpy( rfield.chSpType[rfield.nShape], "READ " );
01112
01113 ReadTable(chFile);
01114
01115
01116 rfield.tslop[rfield.nShape][0] = 1.;
01117 rfield.tslop[rfield.nShape][1] = 0.;
01118
01119
01120 ++rfield.nShape;
01121 return;
01122 }
01123
01124 else if( p.nMatch("TLUSTY") && !p.nMatch("STAR") )
01125 {
01126
01127 fprintf( ioQQQ, " The TABLE TLUSTY command is no longer supported.\n" );
01128 fprintf( ioQQQ, " Please use TABLE STAR TLUSTY instead. See Hazy for details.\n" );
01129 cdEXIT(EXIT_FAILURE);
01130 }
01131
01132 else if( p.nMatch("RUBI") )
01133 {
01134
01135 for( i=0; i < NRUBIN; i++ )
01136 {
01137 rfield.tNu[rfield.nShape][i].set(tnurbn[i]);
01138 rfield.tslop[rfield.nShape][i] = (realnum)log10(fnurbn[i] /tnurbn[i] );
01139 }
01140 rfield.ncont[rfield.nShape] = NRUBIN;
01141 }
01142
01143
01144
01145 else if( p.nMatch("STAR") )
01146 {
01147 char chMetalicity[5] = "", chODFNew[8], chVaryFlag[7] = "";
01148 bool lgHCa = false, lgHNi = false;
01149 long nval, ndim=0;
01150 double Tlow = -1., Thigh = -1.;
01151 double val[MDIM];
01152
01153
01154 if( p.nMatchErase("1-DI") )
01155 ndim = 1;
01156 else if( p.nMatchErase("2-DI") )
01157 ndim = 2;
01158 else if( p.nMatchErase("3-DI") )
01159 ndim = 3;
01160 else if( p.nMatchErase("4-DI") )
01161 ndim = 4;
01162
01163 if( ndim != 0 )
01164 {
01165
01166 sprintf(chVaryFlag,"%1ld-DIM",ndim);
01167 }
01168
01169
01170 rfield.lgTimeVary[p.m_nqh] = p.nMatch( "TIME" );
01171
01172 static const char table[][2][10] = {
01173 {"Z+1.0 ", "p10"},
01174 {"Z+0.75", "p075"},
01175 {"Z+0.5 ", "p05"},
01176 {"Z+0.4 ", "p04"},
01177 {"Z+0.3 ", "p03"},
01178 {"Z+0.25", "p025"},
01179 {"Z+0.2 ", "p02"},
01180 {"Z+0.1 ", "p01"},
01181 {"Z+0.0 ", "p00"},
01182 {"Z-0.1 ", "m01"},
01183 {"Z-0.2 ", "m02"},
01184 {"Z-0.25", "m025"},
01185 {"Z-0.3 ", "m03"},
01186 {"Z-0.4 ", "m04"},
01187 {"Z-0.5 ", "m05"},
01188 {"Z-0.7 ", "m07"},
01189 {"Z-0.75", "m075"},
01190 {"Z-1.0 ", "m10"},
01191 {"Z-1.3 ", "m13"},
01192 {"Z-1.5 ", "m15"},
01193 {"Z-1.7 ", "m17"},
01194 {"Z-2.0 ", "m20"},
01195 {"Z-2.5 ", "m25"},
01196 {"Z-3.0 ", "m30"},
01197 {"Z-3.5 ", "m35"},
01198 {"Z-4.0 ", "m40"},
01199 {"Z-4.5 ", "m45"},
01200 {"Z-5.0 ", "m50"},
01201 {"Z-INF ", "m99"}
01202 };
01203
01204 strncpy( chMetalicity, "p00", sizeof(chMetalicity) );
01205 for (i=0; i < (long)(sizeof(table)/sizeof(*table)); ++i)
01206 {
01207 if (p.nMatchErase(table[i][0]))
01208 {
01209 strncpy( chVaryFlag, table[i][0], sizeof(chVaryFlag) );
01210 strncpy( chMetalicity, table[i][1], sizeof(chMetalicity) );
01211 break;
01212 }
01213 }
01214
01215
01216
01217
01218 bool lgHalo = p.nMatch("HALO");
01219 bool lgSolar = ( !lgHalo && strcmp( chMetalicity, "p00" ) == 0 );
01220
01221
01222 bool lgPG1159 = p.nMatchErase("PG1159");
01223
01224 bool lgList = p.nMatch("LIST");
01225
01226 if( p.nMatch("AVAI") )
01227 {
01228 AtmospheresAvail();
01229 cdEXIT(EXIT_SUCCESS);
01230 }
01231
01232
01233 for( nval=0; nval < MDIM; nval++ )
01234 {
01235 val[nval] = p.FFmtRead();
01236 if( p.lgEOL() )
01237 break;
01238 }
01239
01240 if( nval == 0 && !lgList )
01241 {
01242 fprintf( ioQQQ, " The stellar temperature MUST be entered.\n" );
01243 cdEXIT(EXIT_FAILURE);
01244 }
01245
01246
01247
01248
01249
01250
01251
01252 if( ( val[0] <= 10. && !p.nMatch("LINE") ) || p.nMatch("LOG ") )
01253 {
01254 if( val[0] < log10(BIGFLOAT) )
01255 val[0] = pow(10.,val[0]);
01256 else
01257 {
01258 fprintf(ioQQQ," DISASTER the log of the temperature was specified, "
01259 "but the numerical value is too large.\n Sorry.\n\n");
01260 cdEXIT(EXIT_FAILURE );
01261 }
01262 }
01263
01264 if( lgQuoteFound )
01265 {
01266 nstar = GridInterpolate(val,&nval,&ndim,chFile,lgList,&Tlow,&Thigh);
01267 }
01268
01269 else if( p.nMatch("ATLA") )
01270 {
01271
01272
01273
01274
01275 if( p.nMatch("ODFN") )
01276 strncpy( chODFNew, "_odfnew", sizeof(chODFNew) );
01277 else
01278 strncpy( chODFNew, "", sizeof(chODFNew) );
01279
01280
01281 nstar = AtlasInterpolate(val,&nval,&ndim,chMetalicity,chODFNew,lgList,&Tlow,&Thigh);
01282 }
01283
01284 else if( p.nMatch("COST") )
01285 {
01286
01287
01288
01289
01290 if( p.nMatch("INDE") )
01291 {
01292 imode = IM_COSTAR_TEFF_MODID;
01293 if( nval == 1 )
01294 {
01295 val[1] = 1.;
01296 nval++;
01297 }
01298
01299
01300
01301 if( val[1] < 1. || val[1] > 7. )
01302 {
01303 fprintf( ioQQQ, " The second number must be the id sequence number, 1 to 7.\n" );
01304 fprintf( ioQQQ, " reset to 1.\n" );
01305 val[1] = 1.;
01306 }
01307 }
01308 else if( p.nMatch("ZAMS") )
01309 {
01310 imode = IM_COSTAR_MZAMS_AGE;
01311 if( nval == 1 )
01312 {
01313 fprintf( ioQQQ, " A second number, the age of the star, must be present.\n" );
01314 cdEXIT(EXIT_FAILURE);
01315 }
01316 }
01317 else if( p.nMatch(" AGE") )
01318 {
01319 imode = IM_COSTAR_AGE_MZAMS;
01320 if( nval == 1 )
01321 {
01322 fprintf( ioQQQ, " A second number, the ZAMS mass of the star, must be present.\n" );
01323 cdEXIT(EXIT_FAILURE);
01324 }
01325 }
01326 else
01327 {
01328 if( nval == 1 )
01329 {
01330 imode = IM_COSTAR_TEFF_MODID;
01331
01332 val[1] = 1.;
01333 nval++;
01334 }
01335 else
01336 {
01337
01338
01339
01340 imode = IM_COSTAR_TEFF_LOGG;
01341 }
01342 }
01343
01344 if( ! ( lgSolar || lgHalo ) )
01345 {
01346 fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" );
01347 cdEXIT(EXIT_FAILURE);
01348 }
01349
01350 nstar = CoStarInterpolate(val,&nval,&ndim,imode,lgHalo,lgList,&Tlow,&Thigh);
01351 }
01352
01353 else if( p.nMatch("KURU") )
01354 {
01355 nstar = Kurucz79Interpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01356 }
01357
01358 else if( p.nMatch("MIHA") )
01359 {
01360 nstar = MihalasInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01361 }
01362
01363 else if( p.nMatch("RAUC") )
01364 {
01365 if( ! ( lgSolar || lgHalo ) )
01366 {
01367 fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" );
01368 cdEXIT(EXIT_FAILURE);
01369 }
01370
01371
01372
01373
01374 if( p.nMatch("H-CA") || p.nMatch(" OLD") )
01375 {
01376 lgHCa = true;
01377 nstar = RauchInterpolateHCa(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh);
01378 }
01379 else if( p.nMatch("HYDR") )
01380 {
01381 nstar = RauchInterpolateHydr(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01382 }
01383 else if( p.nMatch("HELI") )
01384 {
01385 nstar = RauchInterpolateHelium(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01386 }
01387 else if( p.nMatch("H+HE") )
01388 {
01389 nstar = RauchInterpolateHpHe(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01390 }
01391 else if( lgPG1159 )
01392 {
01393 nstar = RauchInterpolatePG1159(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01394 }
01395 else if( p.nMatch("CO W") )
01396 {
01397 nstar = RauchInterpolateCOWD(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01398 }
01399 else
01400 {
01401 lgHNi = true;
01402 nstar = RauchInterpolateHNi(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh);
01403 }
01404 }
01405
01406 else if( p.nMatch("TLUS") )
01407 {
01408 if( p.nMatch("OBST") )
01409 {
01410
01411
01412 nstar = TlustyInterpolate(val,&nval,&ndim,TL_OBSTAR,chMetalicity,lgList,&Tlow,&Thigh);
01413 }
01414 else if( p.nMatch("BSTA") )
01415 {
01416
01417
01418 nstar = TlustyInterpolate(val,&nval,&ndim,TL_BSTAR,chMetalicity,lgList,&Tlow,&Thigh);
01419 }
01420 else if( p.nMatch("OSTA") )
01421 {
01422
01423
01424 nstar = TlustyInterpolate(val,&nval,&ndim,TL_OSTAR,chMetalicity,lgList,&Tlow,&Thigh);
01425 }
01426 else
01427 {
01428 fprintf( ioQQQ, " There must be a third key on TABLE STAR TLUSTY;" );
01429 fprintf( ioQQQ, " the options are OBSTAR, BSTAR, OSTAR.\n" );
01430 cdEXIT(EXIT_FAILURE);
01431 }
01432 }
01433
01434 else if( p.nMatch("WERN") )
01435 {
01436
01437
01438 nstar = WernerInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01439 }
01440
01441 else if( p.nMatch("WMBA") )
01442 {
01443
01444
01445 nstar = WMBASICInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh);
01446 }
01447
01448 else
01449 {
01450 fprintf( ioQQQ, " There must be a second key on TABLE STAR;" );
01451 fprintf( ioQQQ, " the options are ATLAS, KURUCZ, MIHALAS, RAUCH, WERNER, and WMBASIC.\n" );
01452 cdEXIT(EXIT_FAILURE);
01453 }
01454
01455
01456 strcpy( rfield.chSpType[rfield.nShape], "VOLK " );
01457
01458 rfield.ncont[rfield.nShape] = nstar;
01459
01460
01461 if( optimize.lgVarOn )
01462 {
01463 optimize.vincr[optimize.nparm] = (realnum)0.1;
01464
01465 if( lgQuoteFound )
01466 {
01467
01468 optimize.nvarxt[optimize.nparm] = nval;
01469 sprintf( optimize.chVarFmt[optimize.nparm], "TABLE STAR \"%s\"", chFile );
01470 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG" );
01471 for( i=1; i < nval; i++ )
01472 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
01473 }
01474
01475 if( p.nMatch("ATLA") )
01476 {
01477
01478 optimize.nvarxt[optimize.nparm] = ndim;
01479 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR ATLAS " );
01480 strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag );
01481 if( p.nMatch("ODFN") )
01482 strcat( optimize.chVarFmt[optimize.nparm], " ODFNEW" );
01483
01484 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG %f" );
01485
01486 if( ndim == 3 )
01487 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
01488
01489 }
01490
01491 else if( p.nMatch("COST") )
01492 {
01493
01494 optimize.nvarxt[optimize.nparm] = 2;
01495
01496 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR COSTAR " );
01497 if( lgHalo )
01498 strcat( optimize.chVarFmt[optimize.nparm], "HALO " );
01499
01500 if( imode == IM_COSTAR_TEFF_MODID )
01501 {
01502 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG , INDEX %f" );
01503 }
01504 else if( imode == IM_COSTAR_TEFF_LOGG )
01505 {
01506 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG %f" );
01507 }
01508 else if( imode == IM_COSTAR_MZAMS_AGE )
01509 {
01510 strcat( optimize.chVarFmt[optimize.nparm], "ZAMS %f LOG %f" );
01511 }
01512 else if( imode == IM_COSTAR_AGE_MZAMS )
01513 {
01514 strcat( optimize.chVarFmt[optimize.nparm], "AGE %f LOG %f" );
01515 optimize.vincr[optimize.nparm] = (realnum)0.5;
01516 }
01517 }
01518
01519 else if( p.nMatch("KURU") )
01520 {
01521
01522 optimize.nvarxt[optimize.nparm] = 1;
01523 strcpy( optimize.chVarFmt[optimize.nparm],
01524 "TABLE STAR KURUCZ %f LOG" );
01525 }
01526
01527 else if( p.nMatch("MIHA") )
01528 {
01529
01530 optimize.nvarxt[optimize.nparm] = 1;
01531 strcpy( optimize.chVarFmt[optimize.nparm],
01532 "TABLE STAR MIHALAS %f LOG" );
01533 }
01534
01535 else if( p.nMatch("RAUC") )
01536 {
01537
01538 optimize.nvarxt[optimize.nparm] = ndim;
01539
01540 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR RAUCH " );
01541
01542 if( p.nMatch("HYDR") )
01543 strcat( optimize.chVarFmt[optimize.nparm], "HYDR " );
01544 else if( p.nMatch("HELI") )
01545 strcat( optimize.chVarFmt[optimize.nparm], "HELIUM " );
01546 else if( p.nMatch("H+HE") )
01547 strcat( optimize.chVarFmt[optimize.nparm], "H+HE " );
01548 else if( lgPG1159 )
01549 strcat( optimize.chVarFmt[optimize.nparm], "PG1159 " );
01550 else if( p.nMatch("CO W") )
01551 strcat( optimize.chVarFmt[optimize.nparm], "CO WD " );
01552 else if( lgHCa )
01553 strcat( optimize.chVarFmt[optimize.nparm], "H-CA " );
01554
01555 if( ( lgHCa || lgHNi ) && ndim == 2 )
01556 {
01557 if( lgHalo )
01558 strcat( optimize.chVarFmt[optimize.nparm], "HALO " );
01559 else
01560 strcat( optimize.chVarFmt[optimize.nparm], "SOLAR " );
01561 }
01562
01563 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG %f" );
01564
01565 if( ndim == 3 )
01566 {
01567 if( p.nMatch("H+HE") )
01568 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
01569 else
01570 strcat( optimize.chVarFmt[optimize.nparm], " %f 3-DIM" );
01571 }
01572 }
01573
01574 if( p.nMatch("TLUS") )
01575 {
01576
01577 optimize.nvarxt[optimize.nparm] = ndim;
01578 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR TLUSTY " );
01579 if( p.nMatch("OBST") )
01580 strcat( optimize.chVarFmt[optimize.nparm], "OBSTAR " );
01581 else if( p.nMatch("BSTA") )
01582 strcat( optimize.chVarFmt[optimize.nparm], "BSTAR " );
01583 else if( p.nMatch("OSTA") )
01584 strcat( optimize.chVarFmt[optimize.nparm], "OSTAR " );
01585 else
01586 TotalInsanity();
01587 strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag );
01588 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG %f" );
01589
01590 if( ndim == 3 )
01591 strcat( optimize.chVarFmt[optimize.nparm], " %f" );
01592 }
01593
01594 else if( p.nMatch("WERN") )
01595 {
01596
01597 optimize.nvarxt[optimize.nparm] = 2;
01598 strcpy( optimize.chVarFmt[optimize.nparm],
01599 "TABLE STAR WERNER %f LOG %f" );
01600 }
01601
01602 else if( p.nMatch("WMBA") )
01603 {
01604
01605 optimize.nvarxt[optimize.nparm] = 3;
01606 strcpy( optimize.chVarFmt[optimize.nparm],
01607 "TABLE STAR WMBASIC %f LOG %f %f" );
01608 }
01609
01610 if( rfield.lgTimeVary[p.m_nqh] )
01611 strcat( optimize.chVarFmt[optimize.nparm], " TIME" );
01612
01613
01614 optimize.nvfpnt[optimize.nparm] = input.nRead;
01615
01616 ASSERT( nval <= LIMEXT );
01617
01618
01619 optimize.vparm[0][optimize.nparm] = (realnum)log10(val[0]);
01620 for( i=1; i < nval; i++ )
01621 optimize.vparm[i][optimize.nparm] = (realnum)val[i];
01622
01623
01624 optimize.varang[optimize.nparm][0] = (realnum)log10(Tlow);
01625 optimize.varang[optimize.nparm][1] = (realnum)log10(Thigh);
01626
01627
01628 ++optimize.nparm;
01629 }
01630 }
01631
01632 else if( p.nMatch(" XDR") )
01633 {
01634
01635
01636
01637
01638 i = 0;
01639 rfield.tNu[rfield.nShape][i].set(rfield.emm);
01640 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
01641
01642 ++i;
01643 rfield.tNu[rfield.nShape][i].set(1e3 / EVRYD * 0.95);
01644 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
01645
01646 ++i;
01647 rfield.tNu[rfield.nShape][i].set(1e3 / EVRYD );
01648 rfield.tslop[rfield.nShape][i] = (realnum)log10(1.);
01649 ++i;
01650 rfield.tNu[rfield.nShape][i].set(1e5 / EVRYD );
01651 rfield.tslop[rfield.nShape][i] = (realnum)log10(pow(100.,-0.7));
01652 ++i;
01653 rfield.tNu[rfield.nShape][i].set(1e5 / EVRYD * 1.05);
01654 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
01655 ++i;
01656 rfield.tNu[rfield.nShape][i].set(rfield.egamry);
01657 rfield.tslop[rfield.nShape][i] = (realnum)log10(SMALLFLOAT);
01658
01659 rfield.ncont[rfield.nShape] = i+1;
01660 }
01661
01662 else
01663 {
01664 fprintf( ioQQQ, " There MUST be a keyword on this line. The keys are:"
01665 "AGN, AKN120, CRAB, COOL, DRAINE, HM96, HM05, _ISM, LINE, POWERlaw, "
01666 "READ, RUBIN, STAR, and XDR.\n Sorry.\n" );
01667 cdEXIT(EXIT_FAILURE);
01668 }
01669
01670
01671
01672
01673 if( strcmp(rfield.chSpType[rfield.nShape],"VOLK ") == 0 )
01674 {
01675 ++rfield.nShape;
01676 return;
01677 }
01678
01679
01680
01681 if( lgNoContinuum )
01682 {
01683 return;
01684 }
01685
01686
01687
01688 if( rfield.tNu[rfield.nShape][0].Ryd() >= 5. )
01689 {
01690 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
01691 {
01692 rfield.tNu[rfield.nShape][i].set(
01693 pow(10.,rfield.tNu[rfield.nShape][i].Ryd() - 15.51718));
01694 }
01695 }
01696
01697 else if( rfield.tNu[rfield.nShape][0].Ryd() < 0. )
01698 {
01699
01700 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
01701 {
01702 rfield.tNu[rfield.nShape][i].set(
01703 (realnum)pow(10.,rfield.tNu[rfield.nShape][i].Ryd()));
01704 }
01705 }
01706
01707
01708 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
01709 {
01710 rfield.tFluxLog[rfield.nShape][i] = rfield.tslop[rfield.nShape][i];
01711 }
01712
01713 for( i=0; i < rfield.ncont[rfield.nShape]-1; i++ )
01714 {
01715 rfield.tslop[rfield.nShape][i] =
01716 (realnum)((rfield.tslop[rfield.nShape][i+1] -
01717 rfield.tslop[rfield.nShape][i])/
01718 log10(rfield.tNu[rfield.nShape][i+1].Ryd()/
01719 rfield.tNu[rfield.nShape][i].Ryd()));
01720 }
01721
01722 if( rfield.ncont[rfield.nShape] > 1 && (rfield.ncont[rfield.nShape] + 1 < rfield.nupper) )
01723 {
01724
01725 for( i=rfield.ncont[rfield.nShape]; i < rfield.nupper; i++ )
01726 {
01727 rfield.tNu[rfield.nShape][i].set(0.);
01728 }
01729 }
01730
01731 if( trace.lgConBug && trace.lgTrace )
01732 {
01733 fprintf( ioQQQ, " Table for this continuum; TNU, TFAC, TSLOP\n" );
01734 for( i=0; i < rfield.ncont[rfield.nShape]; i++ )
01735 {
01736 fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNu[rfield.nShape][i].Ryd(),
01737 rfield.tFluxLog[rfield.nShape][i], rfield.tslop[rfield.nShape][i] );
01738 }
01739 }
01740
01741
01742
01743
01744 lgHit = false;
01745
01746 fac = -DBL_MAX;
01747
01748
01749 for( i=0; i < rfield.ncont[rfield.nShape] && !lgHit; i++ )
01750 {
01751 if( rfield.tNu[rfield.nShape][i].Ryd() > 1. )
01752 {
01753 fac = rfield.tFluxLog[rfield.nShape][i];
01754 lgHit = true;
01755 }
01756 }
01757
01758 if( rfield.ncont[rfield.nShape] > 1 && lgHit )
01759 {
01760
01761 for( i=0; i < rfield.ncont[rfield.nShape] ; i++ )
01762 {
01763 rfield.tFluxLog[rfield.nShape][i] -= (realnum)fac;
01764 }
01765 }
01766
01767 if( rfield.ncont[rfield.nShape] > 1 )
01768 ++rfield.nShape;
01769 return;
01770 }
01771
01772
01773
01774
01775 STATIC void ZeroContin(void)
01776 {
01777
01778 DEBUG_ENTRY( "ZeroContin()" );
01779
01780
01781
01782 tnudrn[0] = 0.3676;
01783 tnudrn[1] = 0.4144;
01784 tnudrn[2] = 0.4558;
01785 tnudrn[3] = 0.5064;
01786 tnudrn[4] = 0.5698;
01787 tnudrn[5] = 0.6511;
01788 tnudrn[6] = 0.7012;
01789 tnudrn[7] = 0.7597;
01790 tnudrn[8] = 0.8220;
01791 tnudrn[9] = 0.9116;
01792 tnudrn[10] = 0.9120;
01793 tnudrn[11] = 0.9306;
01794 tnudrn[12] = 0.9600;
01795 tnudrn[13] = 0.9806;
01796
01797
01798
01799
01800 tnudrn[14] = 0.9999;
01801
01802
01803
01804
01805 int i;
01806 i= 0;
01807 tsldrn[i] = -17.8063;
01808 ++i;tsldrn[i] = -17.7575;
01809 ++i;tsldrn[i] = -17.7268;
01810 ++i;tsldrn[i] = -17.7036;
01811 ++i;tsldrn[i] = -17.6953;
01812 ++i;tsldrn[i] = -17.7182;
01813 ++i;tsldrn[i] = -17.7524;
01814 ++i;tsldrn[i] = -17.8154;
01815 ++i;tsldrn[i] = -17.9176;
01816 ++i;tsldrn[i] = -18.1675;
01817 ++i;tsldrn[i] = -18.1690;
01818 ++i;tsldrn[i] = -18.2477;
01819 ++i;tsldrn[i] = -18.4075;
01820 ++i;tsldrn[i] = -18.5624;
01821 ++i;tsldrn[i] = -18.7722;
01822
01823
01824
01825
01826 tnuagn[0] = 1e-5;
01827 tnuagn[1] = 9.12e-3;
01828 tnuagn[2] = .206;
01829 tnuagn[3] = 1.743;
01830 tnuagn[4] = 4.13;
01831 tnuagn[5] = 26.84;
01832 tnuagn[6] = 7353.;
01833 tnuagn[7] = 7.4e6;
01834
01835 tslagn[0] = -3.388;
01836 tslagn[1] = 4.0115;
01837 tslagn[2] = 2.6576;
01838 tslagn[3] = 2.194;
01839 tslagn[4] = 1.819;
01840 tslagn[5] = -.6192;
01841 tslagn[6] = -2.326;
01842 tslagn[7] = -7.34;
01843 resetBltin( tnuagn , tslagn , true );
01844
01845
01846
01847
01848 tnucrb[0] = 1.0e-5;
01849 tnucrb[1] = 5.2e-4;
01850 tnucrb[2] = 1.5e-3;
01851 tnucrb[3] = 0.11;
01852 tnucrb[4] = 0.73;
01853 tnucrb[5] = 7.3;
01854 tnucrb[6] = 73.;
01855 tnucrb[7] = 7300.;
01856 tnucrb[8] = 1.5e6;
01857 tnucrb[9] = 7.4e6;
01858
01859 fnucrb[0] = 3.77e-21;
01860 fnucrb[1] = 1.38e-21;
01861 fnucrb[2] = 2.10e-21;
01862 fnucrb[3] = 4.92e-23;
01863 fnucrb[4] = 1.90e-23;
01864 fnucrb[5] = 2.24e-24;
01865 fnucrb[6] = 6.42e-26;
01866 fnucrb[7] = 4.02e-28;
01867 fnucrb[8] = 2.08e-31;
01868 fnucrb[9] = 1.66e-32;
01869 resetBltin( tnucrb , fnucrb , false );
01870
01871
01872
01873
01874 resetBltin( tnurbn , fnurbn , false );
01875
01876
01877 cfle[0] = 0.0000100;
01878 cflf[0] = -0.8046910;
01879 cfle[1] = 0.7354023;
01880 cflf[1] = -0.8046910;
01881 cfle[2] = 1.4708046;
01882 cflf[2] = -0.7436830;
01883 cfle[3] = 2.2062068;
01884 cflf[3] = -0.6818757;
01885 cfle[4] = 2.9416091;
01886 cflf[4] = -0.7168990;
01887 cfle[5] = 3.6770115;
01888 cflf[5] = -0.8068384;
01889 cfle[6] = 4.4124136;
01890 cflf[6] = -0.6722584;
01891 cfle[7] = 5.1478162;
01892 cflf[7] = -0.7626385;
01893 cfle[8] = 5.8832183;
01894 cflf[8] = -1.0396487;
01895 cfle[9] = 6.6186204;
01896 cflf[9] = -0.7972314;
01897 cfle[10] = 7.3540230;
01898 cflf[10] = -0.9883416;
01899 cfle[11] = 14.7080460;
01900 cflf[11] = -1.1675659;
01901 cfle[12] = 22.0620689;
01902 cflf[12] = -1.1985949;
01903 cfle[13] = 29.4160919;
01904 cflf[13] = -1.2263466;
01905 cfle[14] = 36.7701149;
01906 cflf[14] = -1.2918345;
01907 cfle[15] = 44.1241379;
01908 cflf[15] = -1.3510833;
01909 cfle[16] = 51.4781609;
01910 cflf[16] = -1.2715496;
01911 cfle[17] = 58.8321838;
01912 cflf[17] = -1.1098027;
01913 cfle[18] = 66.1862030;
01914 cflf[18] = -1.4315782;
01915 cfle[19] = 73.5402298;
01916 cflf[19] = -1.1327956;
01917 cfle[20] = 147.080459;
01918 cflf[20] = -1.6869649;
01919 cfle[21] = 220.620681;
01920 cflf[21] = -2.0239367;
01921 cfle[22] = 294.160919;
01922 cflf[22] = -2.2130392;
01923 cfle[23] = 367.701141;
01924 cflf[23] = -2.3773901;
01925 cfle[24] = 441.241363;
01926 cflf[24] = -2.5326197;
01927 cfle[25] = 514.7816162;
01928 cflf[25] = -2.5292389;
01929 cfle[26] = 588.3218384;
01930 cflf[26] = -2.8230250;
01931 cfle[27] = 661.8620605;
01932 cflf[27] = -2.9502323;
01933 cfle[28] = 735.4022827;
01934 cflf[28] = -3.0774822;
01935 cfle[29] = 1470.8045654;
01936 cflf[29] = -4.2239799;
01937 cfle[30] = 2206.2067871;
01938 cflf[30] = -5.2547927;
01939 cfle[31] = 2941.6091309;
01940 cflf[31] = -6.2353640;
01941 cfle[32] = 3677.0114746;
01942 cflf[32] = -7.1898708;
01943 cfle[33] = 4412.4135742;
01944 cflf[33] = -8.1292381;
01945 cfle[34] = 5147.8159180;
01946 cflf[34] = -9.0594845;
01947 cfle[35] = 5883.2182617;
01948 cflf[35] = -9.9830370;
01949 cfle[36] = 6618.6206055;
01950 cflf[36] = -10.9028034;
01951 cfle[37] = 7354.0229492;
01952 cflf[37] = -11.8188877;
01953 cfle[38] = 7400.0000000;
01954 cflf[38] = -30.0000000;
01955 cfle[39] = 10000000.0000000;
01956 cflf[39] = -30.0000000;
01957 resetBltin( cfle , cflf , true );
01958
01959
01960
01961 tnuakn[0] = 1e-5;
01962 tnuakn[1] = 1.9e-5;
01963 tnuakn[2] = 3.0e-4;
01964 tnuakn[3] = 2.4e-2;
01965 tnuakn[4] = 0.15;
01966 tnuakn[5] = 0.30;
01967 tnuakn[6] = 0.76;
01968 tnuakn[7] = 2.0;
01969 tnuakn[8] = 76.0;
01970 tnuakn[9] = 760.;
01971 tnuakn[10] = 7.4e6;
01972 fnuakn[0] = 1.5e-16;
01973 fnuakn[1] = 1.6e-16;
01974 fnuakn[2] = 1.4e-13;
01975 fnuakn[3] = 8.0e-15;
01976 fnuakn[4] = 1.6e-15;
01977 fnuakn[5] = 1.8e-15;
01978 fnuakn[6] = 7.1e-16;
01979 fnuakn[7] = 7.9e-17;
01980 fnuakn[8] = 1.1e-18;
01981 fnuakn[9] = 7.1e-20;
01982 fnuakn[10] = 1.3e-24;
01983 resetBltin( tnuakn , fnuakn , false );
01984
01985
01986
01987
01988 tnuism[0] = 6.00;
01989
01990 tnuism[1] = 10.72;
01991 tnuism[2] = 11.00;
01992 tnuism[3] = 11.23;
01993 tnuism[4] = 11.47;
01994 tnuism[5] = 11.55;
01995 tnuism[6] = 11.85;
01996 tnuism[7] = 12.26;
01997 tnuism[8] = 12.54;
01998 tnuism[9] = 12.71;
01999 tnuism[10] = 13.10;
02000 tnuism[11] = 13.64;
02001 tnuism[12] = 14.14;
02002 tnuism[13] = 14.38;
02003 tnuism[14] = 14.63;
02004 tnuism[15] = 14.93;
02005 tnuism[16] = 15.08;
02006 tnuism[17] = 15.36;
02007 tnuism[18] = 15.43;
02008 tnuism[19] = 16.25;
02009 tnuism[20] = 17.09;
02010 tnuism[21] = 18.00;
02011 tnuism[22] = 23.00;
02012
02013 fnuism[0] = -16.708;
02014
02015 fnuism[1] = -2.96;
02016 fnuism[2] = -2.47;
02017 fnuism[3] = -2.09;
02018 fnuism[4] = -2.11;
02019 fnuism[5] = -2.34;
02020 fnuism[6] = -3.66;
02021 fnuism[7] = -2.72;
02022 fnuism[8] = -2.45;
02023 fnuism[9] = -2.57;
02024 fnuism[10] = -3.85;
02025 fnuism[11] = -3.34;
02026 fnuism[12] = -2.30;
02027 fnuism[13] = -1.79;
02028 fnuism[14] = -1.79;
02029 fnuism[15] = -2.34;
02030 fnuism[16] = -2.72;
02031 fnuism[17] = -2.55;
02032 fnuism[18] = -2.62;
02033 fnuism[19] = -5.68;
02034 fnuism[20] = -6.45;
02035 fnuism[21] = -6.30;
02036 fnuism[22] = -11.3;
02037
02038 return;
02039 }
02040
02041
02042 int lines_table()
02043 {
02044 DEBUG_ENTRY( "lines_table()" );
02045
02046 if( chLINE_LIST.empty() )
02047 return 0;
02048
02049 vector<char*> chLabel;
02050 vector<realnum> wl;
02051
02052 long nLINE_TABLE = cdGetLineList( chLINE_LIST.c_str(), chLabel, wl );
02053
02054
02055 if( nLINE_TABLE == 0 )
02056 return 0;
02057
02058 fprintf( ioQQQ , "lines_table checking lines within data table %s\n", chLINE_LIST.c_str() );
02059 long miss = 0;
02060
02061
02062 for( long n=0; n < nLINE_TABLE; ++n )
02063 {
02064 double relative , absolute;
02065 if( (cdLine( chLabel[n], wl[n] , &relative , &absolute )) <= 0 )
02066 {
02067 ++miss;
02068 fprintf(ioQQQ,"lines_table in parse_table.cpp did not find line %4s ",chLabel[n]);
02069 prt_wl(ioQQQ,wl[n]);
02070 fprintf(ioQQQ,"\n");
02071 }
02072 }
02073 if( miss )
02074 {
02075
02076 fprintf( ioQQQ , " BOTCHED MONITORS!!! Botched Monitors!!! lines_table could not find a total of %li lines\n\n", miss );
02077 }
02078 else
02079 {
02080 fprintf( ioQQQ , "lines_table found all lines\n\n" );
02081 }
02082
02083
02084 for( unsigned j=0; j < chLabel.size(); ++j )
02085 delete[] chLabel[j];
02086 chLabel.clear();
02087
02088 return miss;
02089 }
02090
02091
02092 STATIC void ReadTable(const char *fnam)
02093 {
02094 char chLine[INPUT_LINE_LENGTH];
02095 long int i;
02096 FILE *io;
02097
02098 DEBUG_ENTRY( "ReadTable()" );
02099
02100
02101 for( i=0; i < NCELL; i++ )
02102 {
02103 rfield.tFluxLog[rfield.nShape][i] = -70.;
02104 }
02105
02106 io = open_data( fnam, "r", AS_LOCAL_ONLY );
02107
02108 string unit;
02109 char *last;
02110
02111
02112 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
02113 {
02114 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02115 goto error;
02116 }
02117
02118 unit = "Ryd";
02119 last = strchr_s(chLine,'\t');
02120 if (last)
02121 {
02122 *last = '\0';
02123 char *first = strrchr(chLine,'/');
02124 if (first)
02125 {
02126 unit = string(first+1);
02127 }
02128 *last = '\t';
02129 };
02130
02131
02132 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
02133 {
02134 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02135 goto error;
02136 }
02137
02138
02139 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
02140 {
02141 long version;
02142 sscanf( chLine, "%ld", &version );
02143 if( version != VERSION_TRNCON )
02144 {
02145 fprintf( ioQQQ,
02146 " the input continuum file has the wrong version number, I expected %li and found %li.\n",
02147 VERSION_TRNCON, version);
02148 goto error;
02149 }
02150 }
02151 else
02152 {
02153 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02154 goto error;
02155 }
02156
02157 char md5sum[NMD5];
02158 long nflux;
02159 double mesh_lo, mesh_hi;
02160 union {
02161 double x;
02162 uint32 i[2];
02163 } u;
02164
02165
02166 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
02167 {
02168 strncpy( md5sum, chLine, NMD5 );
02169 }
02170 else
02171 {
02172 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02173 goto error;
02174 }
02175
02176
02177 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
02178 {
02179 if( cpu.i().big_endian() )
02180 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
02181 else
02182 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
02183 mesh_lo = u.x;
02184 }
02185 else
02186 {
02187 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02188 goto error;
02189 }
02190
02191
02192 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
02193 {
02194 if( cpu.i().big_endian() )
02195 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
02196 else
02197 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
02198 mesh_hi = u.x;
02199 }
02200 else
02201 {
02202 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02203 goto error;
02204 }
02205
02206 if( strncmp( md5sum, continuum.mesh_md5sum.c_str(), NMD5 ) != 0 ||
02207 !fp_equal( mesh_lo, double(rfield.emm) ) ||
02208 !fp_equal( mesh_hi, double(rfield.egamry) ) )
02209 {
02210 fprintf( ioQQQ, " the input continuum file has an incompatible energy grid.\n" );
02211 goto error;
02212 }
02213
02214
02215 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
02216 {
02217 if( cpu.i().big_endian() )
02218 sscanf( chLine, "%x %x", &u.i[0], &u.i[1] );
02219 else
02220 sscanf( chLine, "%x %x", &u.i[1], &u.i[0] );
02221 rfield.RSFCheck[rfield.nShape] = u.x;
02222 }
02223 else
02224 {
02225 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02226 goto error;
02227 }
02228
02229
02230 if( read_whole_line( chLine , (int)sizeof(chLine) , io )!=NULL )
02231 {
02232 sscanf( chLine, "%ld", &nflux );
02233 }
02234 else
02235 {
02236 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02237 goto error;
02238 }
02239
02240
02241 if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
02242 {
02243 fprintf( ioQQQ, " the input continuum file has been truncated.\n" );
02244 goto error;
02245 }
02246
02247
02248 i = 0;
02249
02250 while( (read_whole_line(chLine, (int)sizeof(chLine),io)!=NULL) && (i<NCELL) )
02251 {
02252 double help[2];
02253 sscanf( chLine, "%lf%lf ", &help[0], &help[1] );
02254 rfield.tNu[rfield.nShape][i].set(help[0],unit.c_str());
02255
02256 if (help[1] > 0.0)
02257 {
02258 rfield.tFluxLog[rfield.nShape][i] = (realnum) log10(help[1]/
02259 rfield.tNu[rfield.nShape][i].Ryd());
02260 }
02261 ++i;
02262 }
02263
02264 rfield.ncont[rfield.nShape] = i;
02265
02266
02267 if( rfield.ncont[rfield.nShape] != nflux )
02268 {
02269 fprintf( ioQQQ, " the input continuum file has the wrong number of points: %ld\n",
02270 rfield.ncont[rfield.nShape] );
02271 goto error;
02272 }
02273
02274 fclose(io);
02275 return;
02276
02277 error:
02278 fprintf( ioQQQ, " please recreate this file using the SAVE TRANSMITTED CONTINUUM command.\n" );
02279 cdEXIT(EXIT_FAILURE);
02280 }
02281
02282 #if defined(__HP_aCC)
02283 #pragma OPTIMIZE OFF
02284 #pragma OPTIMIZE ON
02285 #endif