/*GetWerner rebin Werner stellar atmospheres to match cloudy energy grid */ /*DoWernerInterp get a single Werner PN atmosphere, by K Volk */ /*DoWerner read in and interpolate on Werner grid of PN atmospheres, by K Volk */ #include #include "cddefines.h" #include "rfield.h" #include "incidspec.h" #include "path.h" #include "called.h" #include "showme.h" #include "werner.h" /* number of models in the stream */ # ifdef NMODS #undef NMODS # endif #define NMODS 20 /*DoWernerInterp get a single Werner PN atmosphere, by K Volk */ void DoWernerInterp(double temp, double alogg, long int *ierr, FILE * ioIN ); /*RebinWer rebin Werner atmospheres onto Cloudy grid, by Kevin Volk */ void RebinWer( /* the number of points in the incident continuum*/ long n2 , /* the energy grid for the stellar continuum */ float *StarEner , /* the predicted continuum flux on the stellar continuum scale */ float *StarFlux , /* what we want, the continuum on the cloudy grid */ float *CloudyFlux ); /*GetWerner rebin Werner stellar atmospheres to match cloudy energy grid */ void GetWerner(void) { char chLine[100]; long int i, imod; /* these will be malloced into large work arrays*/ float *StarEner , *StarFlux , *CloudyFlux; /* we will write the binary results to this file*/ FILE *ioOUT , /* will get the ascii from this file */ * ioIN ; # ifdef DEBUG_FUN fputs( "<+>GetWerner()\n", debug_fp ); # endif fprintf( ioQQQ, " GetWerner on the job.\n" ); /* malloc some workspace */ StarEner = (float *)malloc( sizeof(float)*NCELL ); if( StarEner == NULL ) { printf( " not enough memory to allocate StarEner in GetWerner\n" ); puts( "[Stop in GetWerner]" ); exit(1); } StarFlux = (float *)malloc( sizeof(float)*NCELL ); if( StarFlux == NULL ) { printf( " not enough memory to allocate StarFlux in GetWerner\n" ); puts( "[Stop in GetWerner]" ); exit(1); } CloudyFlux = (float *)malloc( sizeof(float)*NCELL ); if( CloudyFlux == NULL ) { printf( " not enough memory to allocate CloudyFlux in GetWerner\n" ); puts( "[Stop in GetWerner]" ); exit(1); } /* This is a program to re-bin the Werner stellar model spectra to match the * CLOUDY grid. For short wavelengths I will use a power law extrapolation * of the model values (which should be falling rapidly) if needed. At long * wavelengths I will assume Rayleigh-Jeans from the last stellar model point * to extrapolate to 1 cm wavelength. */ /* This version uses power-law interpolation between the points of the stellar * model. */ /* If ncell is changed the format statement below must also be changed. * If the size of the CLOUDY frequency grid becomes larger than ncell the * arrays below would have to be increased in size. */ /* This file holds the energy bin values (center energy and bin width, both in * Rydberg units) as written out by CLOUDY. The file has two header lines. * * For putting this in CLOUDY the arrrays "anu" and "widflx" would have to be * passed into the routine. */ /* This "werner.ascii" file is a straight ascii dump of the Klaus Werner * stellar model files which he gave to me in 1992. The first set of values * is the 80000 K log(g)=5.0 model fluxes, followed by the wavelength grid, * followed by the other models in order of increasing log(g) and temperature. */ if( (ioIN = fopen( "werner.ascii", "r" ) ) == NULL ) { fprintf( ioQQQ, " GetWerner fails opening werner.ascii\n" ); fprintf( ioQQQ, " Where is this file?\n" ); puts( "[Stop in getwerner]" ); exit(1); } if( (ioOUT = fopen( "werner.mod", "wb" ) ) == NULL ) { fprintf( ioQQQ, " GetWerner fails creating werner.mod\n" ); fprintf( ioQQQ, " This is impossible??\n" ); puts( "[Stop in getwerner]" ); exit(1); } fprintf( ioQQQ, " GetWerner creating werner.mod\n" ); /* first dump the Cloudy energy scale in the binary file. This will * be used for sanity checks when the file is read*/ if( abs( fwrite( rfield.AnuOrg, 1,sizeof(rfield.AnuOrg),ioOUT ) - sizeof(rfield.AnuOrg) ) !=0 ) { fprintf( ioQQQ, " problem1 trying to write werner.mod\n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(rfield.AnuOrg) ); puts( "[Stop in GetWerner]" ); exit(1); } /* ======================================================================== */ /* get energy grid, */ i = 0; while( i < 514 ) { if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL ) { fprintf( ioQQQ, " GetWerner fails reading energy grid1.\n" ); puts( "[Stop in GetData]" ); exit(-1); } sscanf( chLine , "%f %f %f %f %f" , &StarEner[i],&StarEner[i+1],&StarEner[i+2],&StarEner[i+3],&StarEner[i+4] ); /* increment counter, five numbers per line*/ i += 5; } for( i=1; i<514; ++i) { /* following is due to there being fluxes above and below certain * wavelengths where the opacity changes * (i.e. the Lyman and Balmer limits for example) which are * assigned the same wavelength in the Klaus Werner files. */ if( StarEner[i] == StarEner[i-1] ) { StarEner[i-1] *= 0.999f; StarEner[i] *= 1.001f; } } /* The Klaus Werner model wavelength grid was converted to Rydberg units * when I originally created the direct access file, but here they are * converted back to Hz. */ /* ======================================================================== */ for( imod=1; imod < 21; imod++ ) { /* now get the models */ i = 0; while( i < 514 ) { if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL ) { fprintf( ioQQQ, " GetWerner fails reading star file1.\n" ); puts( "[Stop in GetData]" ); exit(-1); } sscanf( chLine , "%f %f %f %f %f" , &StarFlux[i],&StarFlux[i+1],&StarFlux[i+2], &StarFlux[i+3],&StarFlux[i+4] ); /* increment counter, five numbers per line*/ i += 5; } /* continuum flux was log, convert to linear, * and make sure no doubled energy points */ for( i=0; i < 514; i++ ) { StarFlux[i] = (float)pow(10.,StarFlux[i]); } /* this will do the heavy lifting, and define arrays used below * 514 is the number of energy points*/ RebinWer(514, StarEner, StarFlux, CloudyFlux ); /*{ FILE *ioBUG; ioBUG = fopen("test.txt","w"); for( i=0; iGetWerner()\n", debug_fp ); # endif return; } /*====================================================================== */ /*RebinWer rebin Werner atmospheres onto Cloudy grid, by Kevin Volk */ void RebinWer( /* the number of points in the incident continuum*/ long n2 , /* the energy grid for the stellar continuum */ float *StarEner , /* the predicted continuum flux on the stellar continuum scale */ float *StarFlux , /* what we want, the continuum on the cloudy grid */ float *CloudyFlux ) { long int i, ipHi, ipLo, j; double EdgeHigh, EdgeLow, power, slope2, sum, v1, /*v2, */ val, x1, x2; # ifdef DEBUG_FUN fputs( "<+>RebinWer()\n", debug_fp ); # endif /* this is the Rayleigh-Jeans slope in F_nu, for low frequency extrapolation */ slope2 = 2.0e00; for( j=0; j < NCELL; j++ ) { /* EdgeHigh is upper bound of this continnum cell */ EdgeHigh = rfield.anu[j] + 0.5e00*rfield.widflx[j]; /* EdgeLow is lower bound of this continnum cell */ EdgeLow = rfield.anu[j] - 0.5e00*rfield.widflx[j]; if( EdgeHigh < StarEner[0] ) { /* this is case where cloudy's continuum is below stellar continuum, * so we do Rayleigh Jeans extrapolation */ CloudyFlux[j] = (float)(StarFlux[0]*(pow(rfield.anu[j]/ StarEner[0],slope2))); } else { if( EdgeLow > StarEner[n2-1] ) { /* case where cloudy continuum above highest stellar point */ CloudyFlux[j] = 0.0e00; } else { /* now go through stellar continuum to find bins corresponding to * this cloudy cell, stellar continuum defined through n2 cells */ ipLo = 0; ipHi = 0; for( i=0; i < n2; i++ ) { if( StarEner[i+1] > EdgeLow && StarEner[i] <= EdgeLow ) { ipLo = i; break; } } for( i=MAX2(0,ipLo-1); i < n2; i++ ) { if( StarEner[i+1] > EdgeHigh && StarEner[i] <= EdgeHigh ) { ipHi = i; break; } } /* Do the case where the grid point and its edges are between two of the * stellar model points: add segments with power-law interpolation up to * do the averaging.*/ /* When EdgeLow < anu3(1) ipLo is equal to -1 and ipHi is equal to 1, so do it * as if ipHi = ipLo to avoid going outside array bounds. */ if( ipHi == ipLo || ipLo>ipHi) { power = log(StarFlux[ipLo+1]/StarFlux[ipLo])/log(StarEner[ipLo+1]/ StarEner[ipLo]); val = power*log(rfield.anu[j]/StarEner[ipLo]); CloudyFlux[j] = (float)(StarFlux[ipLo]*exp(val)); } else { sum = 0.0e00; /* ipHi points to stellar point at high end of cloudy continuum cell, * ipLo points to low end */ for( i=ipLo; i <= ipHi; i++ ) { power = log(StarFlux[i+1]/StarFlux[i])/ log(StarEner[i+1]/StarEner[i]); if( i == ipLo ) { x1 = EdgeLow; x2 = StarEner[i+1]; v1 = StarFlux[i]*exp(power*log(x1/StarEner[i])); /*v2 = StarFlux[i+1];*/ } else if( i == ipHi ) { x2 = EdgeHigh; x1 = StarEner[i]; /*v2 = StarFlux[i]*exp(power*log(x1/StarEner[i]));*/ v1 = StarFlux[i]; } /*if( i > ipLo && i < ipHi )*/ else { x1 = StarEner[i]; x2 = StarEner[i+1]; v1 = StarFlux[i]; /*v2 = StarFlux[i+1];*/ } if( fabs(power+1.0) < 0.001 ) { val = x1*v1*log(x2/x1); } else { val = pow(x2/x1,power + 1.e00) - 1.e00; val = val*x1*v1/(power + 1.e00); } sum += val; } CloudyFlux[j] = (float)(sum/rfield.widflx[j]); } } } } # ifdef DEBUG_FUN fputs( " <->RebinWer()\n", debug_fp ); # endif return; } /*DoWerner read in and interpolate on Werner grid of PN atmospheres, by K Volk */ void DoWerner(long int *nstar, double temp, double alogg) { /* will be used for reading in werner.mod */ FILE *ioIN ; char chLine[81]; int lgAgain; long int ierr, j, jmod; double dev1, dev2; float wl[NCELL]; static float teff[NMODS]={80000.,80000.,80000.,80000.,100000., 100000.,100000.,100000.,120000.,120000.,120000.,140000.,140000., 140000.,160000.,160000.,180000.,180000.,200000.,200000.}; static float xlogg[NMODS]={5.,6.,7.,8.,5.,6.,7.,8.,6.,7.,8.,6., 7.,8.,7.,8.,7.,8.,7.,8.}; # ifdef DEBUG_FUN fputs( "<+>DoWerner()\n", debug_fp ); # endif /* structure of werner.mod * * all recordes are floats NCELL long, * the first record contaings an image of the cloudy energy * grid, to use for sanity checks * the remaining records are rebinned stellar continua * there are a total of 20 of these */ /* parameter (ndata=865,nmods=20) */ /* SUBROUTINE "WERNE" WAS ADDED (28 DEC 1992) TO READ FROM THE SET OF * HOT WHITE DWARF MODEL ATMOSPHERES FROM KLAUS WERNER AT KEIL. THE * VALUES ARE READ IN (ENERGY IN RYDBERG UNITS, LOG10(F_NU) [F_NU IN CGS * UNITS]) FOR ANY OF 20 MODELS. T_EFF VALUES OF 80, 100, 120, 140, 160, * 180 AND 200 THOUSAND K AND LOG(G) VALUES OF 5.0, 6.0, 7.0, AND 8.0 ARE * AVAILABLE. THIS VERSION DOES INTERPOLATION. IF A MODEL HAS A T_EFF * AND A G VALUE WITHIN 1% OF THE VALUES SPECIFIED IN THE PARAMETERS, THAT * MODEL IS CHOSEN FOR USE. EACH MODEL HAS 514 POINTS, THE LONGEST WAVELENGTH * POINT BEING AN EXTRAPOLATION FROM THE NEXT TWO POINTS. THE FILE IS A * FORMATTED DIRECT ACCESS FILE WITH RECORD LENGTH 7*865=6055 BYTES, FORMAT * 865I7=3*250i7,115i7. * * FILE USED: WERNER.MOD (CONTAINS THE WAVELENGTH ARRAY, THEN THE MODELS) * * TESTS INDICATE THAT THE INTERPOLATION AND RESCALING PRESERVE THE SET * T_EFF VALUE (MEASURED BY THE SURFACE FLUX) TO WITHIN ABOUT 100 K FOR * ANY TEMPERATURE IN THE 80,000 TO 200,000 K RANGE. RESCALING BY FACTORS * OF ORDER 1.01 TO 1.1 ARE NEEDED TO ACHIEVE THIS, SINCE THE INTERPOLATION BY * ITSELF CAN PRODUCE SOME DEVIATIONS. THE LARGEST RESCALING IS NEEDED AT * "LOW" T_EFF AND LOG(G) VALUES. FOR HIGH TEMPERATURES THE SCALING IS SELDOM * LARGER THAN 1.01. * * THESE RESCALING RESULTS ARE FOR THE FULL MODELS, NOT THE REBINNED VERSIONS. * THE NUMBER OF GRID POINTS IS NOT CHANGED TOO MUCH, SO THE RESULTS SHOULD * NOT BE AFFECTED TOO MUCH. * * * nmods is the number of accessable models in the werner model set. * data set 1 is the wavelength grid, so there are nmods+1 data sets * in the main file. each set has 865 re-binned values, increasing in * frequency. the original models had 514 values. * */ /* true is special path set, false if file in same dir */ if( lgDataPathSet ) { /* path is set, generate full path name with file */ strcpy( chLine , chDataPath ); strcat( chLine, "werner.mod" ); } else { /* path not set, look here */ strcpy( chLine , "werner.mod" ); } if( (ioIN = fopen( chLine , "rb" )) == NULL ) { /* something went wrong */ fprintf( ioQQQ, "ERROR: The Klaus Werner stellar atmosphere file was not found.\nSorry.\n" ); puts( "[Stop in DoWerner]" ); exit(1); } /* read in the saved cloudy energy scale so we can confirm this is a good image */ if( abs( fread( wl, 1, sizeof(wl),ioIN ) - sizeof(wl) ) !=0 ) { fprintf( ioQQQ, " problem trying to read werner.mod wavelengths \n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(wl) ); puts( "[Stop in DoWerner]" ); exit(1); } /* TEMP and ALOGG are desired temp and log g */ /* this may (or may not) catch bad input parameters */ if( temp < 80000. || temp > 200000. ) { fprintf( ioQQQ, " Entered temperature is outside range of grid.\n" ); puts( "[Stop in DoWerner]" ); exit(1); } if( alogg < 5. || alogg > 8. ) { fprintf( ioQQQ, " Entered gravity is outside range of grid.\n" ); puts( "[Stop in DoWerner]" ); exit(1); } lgAgain = TRUE; j = 0; /* -1 will be indication that interpolation is needed */ jmod = -1; while( lgAgain && j < NMODS ) { dev1 = fabs(temp/teff[j]-1.); dev2 = fabs(alogg-xlogg[j]); /* this is a test for how close we are to the tabulated models */ if( dev1 <= 0.01 && dev2 <= 0.00414 ) { /* jmod is pointer to record within werner.mod file */ /* the +1 below is to offset the fseek to one beyond * the hit in the array - this is to skip the set of energies, * which are store at the first file but not included in * the list of models */ jmod = j; lgAgain = FALSE; } ++j; } /* if jmod still -1 then we need to interpolate */ if( jmod == -1 ) { /* was called getmod but name changed to following */ DoWernerInterp( temp , alogg , &ierr, ioIN); if( ierr > 0 ) { fprintf( ioQQQ, " Error occurred in getmod.\n" ); puts( "[Stop in DoWerner]" ); exit(1); } } else { /* this branch, we hit a proper model, use exactly that one */ if( fseek(ioIN, (jmod+1)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", jmod ); puts( "[Stop in DoWerner]" ); exit(1); } if( abs( fread( rfield.flux, 1 , sizeof(rfield.flux), ioIN ) - sizeof(rfield.flux)) != 0 ) { fprintf( ioQQQ, " Error trying to read Werner atmosphere%4ld\n", jmod ); puts( "[Stop in DoWerner]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< K.Volk s K. Werner model %2ld read in. T_eff = %8.1f LOG(g) = %6.3f>>> *\n", jmod+1 , teff[jmod], xlogg[jmod] ); } } for( j=0; j < NCELL ; j++ ) { IncidSpec.tslop[j][IncidSpec.nspec-1] = (float)rfield.flux[j]; /* we save the wavelength scale so that it can be cross checked agains * one inside cloudy - this checks whether the grid is for the current * version of the wavelength scale */ IncidSpec.tnu[j][IncidSpec.nspec-1] = wl[j]; } *nstar = NCELL; fclose( ioIN ); # ifdef DEBUG_FUN fputs( " <->DoWerner()\n", debug_fp ); # endif return; } /*DoWernerInterp interpolate T and G on the Werner PN atmosphere, by K Volk */ void DoWernerInterp( /* desired stellar temperature */ double temp, /* desired surface gravity */ double alogg, /* this will be zero if ended ok */ long int *ierr, /* the file with the stellar continuum */ FILE * ioIN ) { long int ipGrav, ipGravUp , j, ipTeff, ipTeffUp; /* following two used in binary read, must remain float * or binary file changed with it */ float flux1[NCELL], flux2[NCELL]; double fact1, fact2, fr1, fr2, fr3, fr4; /* offset for fseek to first star file */ const long iOff=-1; /* the computed gravities */ static float gval[4]={5.,6.,7.,8.}; /* the computed temperatures*/ static float tval[7]= {80000.,100000.,120000.,140000.,160000.,180000., 200000.}; /* these are pointers to specific models, organized by grav and teff */ static long modnum[4][7]= { {2,6,10,13,16,18,20}, {3,7,10,13,16,18,20}, {4,8,11,14,16,18,20}, {5,9,12,15,17,19,21} }; # ifdef DEBUG_FUN fputs( "<+>DoWernerInterp()\n", debug_fp ); # endif /* this is error indicator, error condition detected if >0 upon exit */ *ierr = 1; /* this subroutine gets stellar atmospheres from the Werner file * and interpolates on them, returning the flux in FLUX as F_nu */ /* pointer to temp will range from 0 to 6, there are seven temps stored */ ipTeff = (long)((temp-80000.)/20000.); ipTeff = MAX2(0, ipTeff ); ipTeff = MIN2( ipTeff , 6); /* upper model, make sure grid not overrun */ ipTeffUp = MIN2( 6, ipTeff + 1 ); /* check whether we can deal with the surface gravity */ /* tabulated alogg is 5 to 8, * pointer will range from 0 to 3 */ ipGrav = (long)(alogg) - 5; ipGrav = MAX2( 0 , ipGrav ); ipGrav = MIN2( 3 , ipGrav ); /* upper model, make sure grid not overrun */ ipGravUp = MIN2( 3, ipGrav + 1 ); /* at this point we should have valid gravity (ipGrav) and * teff (ipTeff) pointers */ /* Read in first model ==================== */ if( fseek(ioIN, (modnum[ipGrav][ipTeff]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", modnum[ipGrav][ipTeff] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( abs( fread( flux1, 1 , sizeof(flux1), ioIN ) - sizeof(flux1)) != 0 ) { fprintf( ioQQQ, " Error trying to read Werner atmosphere%4ld\n", modnum[ipGrav][ipTeff] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } /* say what we got */ if( called.lgTalk ) { fprintf( ioQQQ, " * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", modnum[ipGrav][ipTeff], tval[ipTeff], gval[ipGrav] ); } if( fseek(ioIN, (modnum[ipGravUp][ipTeff]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", modnum[ipGravUp][ipTeff] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } /* Read in second model ==================== */ if( abs( fread( flux2, 1 , sizeof(flux2), ioIN ) - sizeof(flux2)) != 0 ) { fprintf( ioQQQ, " Error2 trying to read Werner atmosphere%4ld\n", modnum[ipGravUp][ipTeff] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< Klaus Werner stellar model %2ld read in." " T_eff = %8.1f log(g) = %6.3f >>> *\n", modnum[ipGravUp][ipTeff] , tval[ipTeff], gval[ipGravUp] ); } fr1 = alogg - gval[ipGrav]; fr2 = 1. - fr1; if( ipTeffUp == ipTeff ) { fr3 = 1.; } else { fr3 = (temp - tval[ipTeff])/(tval[ipTeffUp] - tval[ipTeff]); } fr4 = 1. - fr3; if( fr3 < 0. || fr3 > 1. ) { fprintf( ioQQQ, "fr3 insanity in DoWernerInterp\n" ); ShowMe(); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( fr2 < 0. || fr2 > 1. ) { fprintf( ioQQQ, "fr2 insanity in DoWernerInterp\n" ); ShowMe(); puts( "[Stop in DoWernerInterp]" ); exit(1); } /* save inerpolated continuum in flux array */ for( j=0; j < NCELL; j++ ) { /* this is the interpolation between the two read in so far */ rfield.flux[j] = (float)( fr2*log10(MAX2(1e-37,flux1[j])) + fr1*log10(MAX2(1e-37,flux2[j]))); } /* Read in third model ==================== */ if( fseek(ioIN, (modnum[ipGrav][ipTeffUp]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", modnum[ipGrav][ipTeffUp] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( abs( fread( flux1, 1 , sizeof(flux1), ioIN ) - sizeof(flux1)) != 0 ) { fprintf( ioQQQ, " Error trying to read Werner atmosphere%4ld\n", modnum[ipGrav][ipTeffUp] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", modnum[ipGrav][ipTeffUp] , tval[ipTeffUp], gval[ipGrav] ); } /* Read in fourth model ==================== */ if( fseek(ioIN, (modnum[ipGravUp][ipTeffUp]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", modnum[ipGravUp][ipTeffUp] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( abs( fread( flux2, 1 , sizeof(flux2), ioIN ) - sizeof(flux2)) != 0 ) { fprintf( ioQQQ, " Error3 trying to read Werner atmosphere%4ld\n", modnum[ipGravUp][ipTeffUp] ); puts( "[Stop in DoWernerInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", modnum[ipGravUp][ipTeffUp], tval[ipTeffUp], gval[ipGravUp] ); } /* now do the actual interpolation */ for( j=0; j < NCELL; j++ ) { flux1[j] = (float)( fr2*log10(MAX2(1e-37,flux1[j])) + fr1*log10(MAX2(1e-37,flux2[j]))); } /* now, allow for the total flux being significantly different between the * two set of fluxes since the t_eff values are significantly different. * scale things by the 4th power of t_eff */ fact1 = log10(temp/tval[ipTeff])*4.; fact2 = log10(temp/tval[ipTeffUp])*4.; fact1 = 0.; fact2 = 0.; for( j=0; j < NCELL; j++ ) { rfield.flux[j] = (float)pow(10.,(fr4*(rfield.flux[j] + fact1) + fr3*(flux1[j] + fact2))); if( rfield.flux[j]<1e-36 ) rfield.flux[j] = 0.; } /* say that things turned out ok */ *ierr = 0; # ifdef DEBUG_FUN fputs( " <->DoWernerInterp()\n", debug_fp ); # endif return; }