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