/*GetAtlas rebin Kurucz stellar models to match energy grid of code */ #include #include "cddefines.h" #include "rfield.h" #include "incidspec.h" #include "path.h" #include "called.h" #include "showme.h" #include "atlas.h" #define NATLAS 412 #define NVALS 1221 /*DoAtlasInterp get one of the Atlas model atmospheres, coded by K. Volk */ void DoAtlasInterp(double temp, double alogg, long int jval[][11], float tval[], long int *ierr, float teff[], float xlogg[], FILE *ioIN ); /*RebinAtlas rebin the atlas continuum grid onto the cloudy grid */ void RebinAtlas( /* 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); void GetAtlas(void) { char chLine[100] , chC7[8],chC12[13]; long int i, imod, j; /* these will be malloced into large work arrays*/ float *StarEner , *StarFlux , *CloudyFlux , *scratch; float c, rydb, con1; /* we will write the binary results to this file*/ FILE *ioOUT , /* will get the ascii from this file */ * ioIN , /* this is short list of models in main file */ * ioLIST; # ifdef DEBUG_FUN fputs( "<+>GetAtlas()\n", debug_fp ); # endif /* This is a program to re-bin the Kurucz stellar models spectrum 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.*/ /* nvals is the number of values per Kurucz model. nmods is the total * number of models in the set. The last two are the Vega and Solar models.*/ fprintf( ioQQQ, " GetAtlas on the job.\n" ); /* malloc some workspace */ StarEner = (float *)malloc( sizeof(float)*NCELL ); if( StarEner == NULL ) { printf( " not enough memory to allocate StarEner in GetAtlas\n" ); puts( "[Stop in GetAtlas]" ); exit(1); } StarFlux = (float *)malloc( sizeof(float)*NCELL ); if( StarFlux == NULL ) { printf( " not enough memory to allocate StarFlux in GetAtlas\n" ); puts( "[Stop in GetAtlas]" ); exit(1); } CloudyFlux = (float *)malloc( sizeof(float)*NCELL ); if( CloudyFlux == NULL ) { printf( " not enough memory to allocate CloudyFlux in GetAtlas\n" ); puts( "[Stop in GetAtlas]" ); exit(1); } scratch = (float *)malloc( sizeof(float)*NCELL ); if( scratch == NULL ) { printf( " not enough memory to allocate scratch in GetAtlas\n" ); puts( "[Stop in GetAtlas]" ); exit(1); } /* open main file of stellar atmospheres for reading */ if( (ioIN = fopen( "kurucz.ascii", "r" ) ) == NULL ) { fprintf( ioQQQ, " GetAtlas could not find kurucz.ascii.\n" ); puts( "[Stop in getatlas]" ); exit(1); } fprintf( ioQQQ, " GetAtlas got kurucz.ascii.\n" ); /* open much smaller list of stellar atmospheres */ if( (ioLIST = fopen( "kurucz.list", "r" ) ) == NULL ) { fprintf( ioQQQ, " GetAtlas could not find kurucz.list.\n" ); puts( "[Stop in getatlas]" ); exit(1); } fprintf( ioQQQ, " GetAtlas got kurucz.list.\n" ); /* create the binary output file */ if( (ioOUT = fopen( "atlas.mod", "wb" ) ) == NULL ) { fprintf( ioQQQ, " GetAtlas failed creating kurucz.mod.\n" ); puts( "[Stop in getatlas]" ); exit(1); } fprintf( ioQQQ, " GetAtlas got atlas.mod.\n" ); /* write out the cloudy energy grid for later sanity checks */ if( abs( fwrite( rfield.AnuOrg, 1,sizeof(rfield.AnuOrg),ioOUT ) - sizeof(rfield.AnuOrg) ) !=0 ) { fprintf( ioQQQ, " GetAtlas failed writng anu array.\n" ); puts( "[Stop in getatlas]" ); exit(1); } /* first read wavelenghth grid from start of large data file*/ i = 0; while( i < NVALS ) { if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL ) { fprintf( ioQQQ, " GetAtlas could not read kurucz wavelengths.\n" ); puts( "[Stop in GetAtlas]" ); exit(-1); } /* this is the wavelength grid */ sscanf( chLine , "%f %f %f %f %f" , &scratch[i],&scratch[i+1],&scratch[i+2],&scratch[i+3],&scratch[i+4] ); /* increment counter, five numbers per line*/ i += 5; } fprintf( ioQQQ, "GetAtlas got kurucz wavelengths, low to high is %11.2e to %11.2e\n", scratch[0], scratch[NVALS-1] ); /* the wavelength values scratch are in nm units, convert to freq (Hz) * the energy range is 5.7e-4 to 10 Ryd */ c = 299792458.0f; rydb = 10967758.30f; /* These are the speed of light and H-atom Rydberg constant in MKS units, for * conversion of frequency to rydberg units. The infinite mass * R value is 10973731.534 and the proton to electron mass ratio is * 1836.152701 for this calculation. */ con1 = 1.f/(c*rydb); for( i=0; i < NVALS; i++ ) { /* convert wavelength (in nm) grid to rydberg units */ scratch[i] = (float)(299792458.e09/scratch[i]*con1); } /* flip them over */ for( j=0; j < NVALS; j++ ) { /* StarEner now has correct energy grid, in Rydbergs, in increasing order */ StarEner[j] = scratch[NVALS-j-1]; } /* now get grid of effective temperatures and surface gravities from smaller file */ /* first line is header */ if( fgets( chLine , sizeof(chLine) , ioLIST ) == NULL ) { fprintf( ioQQQ, " GetAtlas could not read empty header.\n" ); puts( "[Stop in GetAtlas]" ); exit(-1); } /* now read in NATLAS pairs of effective temperatures and surface gravities */ for( j=0; j < NATLAS; j++ ) { if( fgets( chLine , sizeof(chLine) , ioLIST ) == NULL ) { fprintf( ioQQQ, " GetAtlas failed readin teff, logg, for model%4ld\n", j ); puts( "[Stop in GetAtlas]" ); exit(-1); } /* parse the temperature and gravity off the line */ sscanf( chLine , "%6li%s%8f %s %f " , &i, chC12, &scratch[j], chC7, &CloudyFlux[j] ); } fprintf( ioQQQ, " Got grid of temperatures and gravities.\n" ); /* dump the temperatures into the binary file */ if( abs( fwrite( scratch , 1, NCELL*sizeof(float ),ioOUT ) - NCELL*sizeof(float ) ) !=0 ) { fprintf( ioQQQ, " GetAtlas failed writing teff.\n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(scratch) ); puts( "[Stop in GetAtlas]" ); exit(1); } /* dump the gravities into the binary file */ if( abs( fwrite( CloudyFlux , 1, NCELL*sizeof(float ),ioOUT ) - NCELL*sizeof(float ) ) !=0 ) { fprintf( ioQQQ, " GetAtlas failed writing grav.\n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(CloudyFlux) ); puts( "[Stop in GetAtlas]" ); exit(1); } /* *** NOTE **** WARNING!!!!! ****** * As things are written here, the CLOUDY grid range covered by the models * must be > 1.0e-05 Rydbergs and < 100000.0 Rydbergs so that the integer * values are not larger than 9,999,999 and are all negative before the abs * function is applied to the value to make the iteff element. 6 digits are * needed in the logarithm value to give 5 digit accuracy when the values are * converted back to energy values in Rydbergs by programs reading the model * file. For the present Werner and Kurucz models these limits are not a * problem, nor are these limits likely to be violated in the future by purely * stellar models....but if the limits are exceeded, things will not work. * * Since the Kurucz models cover up to a temperature of 50000 K and there is * no significant emission at very high energies, the range of energy covered * by the rebinning does not have to extend beyond a few Rydberg. Here 10 Rydb * is used to be consistent with the output from the hotter Werner models.*/ /* write out the cloudy energy grid for later sanity checks */ if( abs( fwrite( rfield.AnuOrg , 1, sizeof(rfield.AnuOrg ),ioOUT ) - sizeof( rfield.AnuOrg ) ) !=0 ) { fprintf( ioQQQ, " GetAtlas failed writing anu array.\n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(rfield.AnuOrg) ); puts( "[Stop in GetAtlas]" ); exit(1); } /* This array will now be used for writing out the flux values in log form, * multiplied by 1.e+05 and converted to integer (this allows sign plus 7 * integer places, storing 5 places in the LOG10 value which gives 4 significant * figure accuracy in the values).*/ fprintf( ioQQQ, "\n\n\n\n" ); fprintf( ioQQQ, " Now would be a good time to get on the Cloudy mailing list, to\n" ); fprintf( ioQQQ, " find out about corrections and improvements to the code.\n" ); fprintf( ioQQQ, " Send a request to cldmail@pa.uky.edu\n" ); fprintf( ioQQQ, "\n\n\n\n" ); fprintf( ioQQQ, " I will print every 10th model so you know I am still alive.\n" ); fprintf( ioQQQ, " (there are 412 models)\n" ); for( imod=0; imod < NATLAS; imod++ ) { i = 0; while( i < NVALS ) { if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL ) { fprintf( ioQQQ, " GetAtlas failed reading star flux.\n" ); puts( "[Stop in GetAtlas]" ); exit(-1); } /* actually scan off the fluxes, 5 per line */ sscanf( chLine , "%f %f %f %f %f" , &scratch[i],&scratch[i+1],&scratch[i+2],&scratch[i+3],&scratch[i+4] ); /* increment the counter */ i += 5; } for( j=0; j < NVALS; j++ ) { StarFlux[j] = scratch[NVALS-1-j]; } /* actually do the rebinning. */ RebinAtlas(NVALS,StarEner, StarFlux, CloudyFlux); /*{ FILE *ioBUG; ioBUG = fopen("test.txt","w"); for( i=0; iGetAtlas()\n", debug_fp ); # endif return; } /*RebinAtlas rebin the atlas continuum grid onto the cloudy grid */ void RebinAtlas( /* 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( "<+>RebinAtlas()\n", debug_fp ); # endif slope2 = 2.0; /* this is the Rayleigh-Jeans slope in F_nu, for low frequency extrapolation*/ /* NCELL is number of cloudy energy cells */ for( j=0; j < NCELL; j++ ) { /* EdgeHigh is upper bound of this continnum cell */ /* lower and upper bounds of this energy cell */ EdgeHigh = rfield.anu[j] + 0.5*rfield.widflx[j]; /* EdgeLow is lower bound of this continnum cell */ EdgeLow = rfield.anu[j] - 0.5*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] ) { CloudyFlux[j] = 0.0; } 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 0 and ipHi is equal to 1, so do it * as if ipHi = ipLo to avoid going outside array bounds. */ if( StarFlux[ipLo]==0. || StarFlux[ipHi+1]==0. ) { CloudyFlux[j] = 0.; } else 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.0; for( i=ipLo; i <= ipHi; i++ ) { /* local continuum slope */ 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.001 ) { val = x1*v1*log(x2/x1); } else { val = pow(x2/x1,power + 1.) - 1.; val = val*x1*v1/(power + 1.); } sum += val; } CloudyFlux[j] = (float)(sum/rfield.widflx[j]); } } } } # ifdef DEBUG_FUN fputs( " <->RebinAtlas()\n", debug_fp ); # endif return; } /*DoAtlas interpolate on atlas model atmospheres, by K Volk */ void DoAtlas(long int *nstar, double temp, double alogg, int lgTrace) { /* will point to binary file */ FILE *ioIN; char chLine[81]; long int i, ierr, j, jj, jmod, jval[61][11], k; double dev1, dev2; float teff[NATLAS], xlogg[NATLAS]; static float tval[61]={3500.,3750.,4000.,4250.,4500.,4750.,5000., 5250.,5500.,5750.,6000.,6250.,6500.,6750.,7000.,7250.,7500., 7750.,8000.,8250.,8500.,8750.,9000.,9250.,9500.,9750.,10000., 10500.,11000.,11500.,12000.,12500.,13000.,14000.,15000.,16000., 17000.,18000.,19000.,20000.,21000.,22000.,23000.,24000.,25000., 26000.,27000.,28000.,29000.,30000.,31000.,32000.,33000.,34000., 35000.,37500.,40000.,42500.,45000.,47500.,50000.}; # ifdef DEBUG_FUN fputs( "<+>DoAtlas()\n", debug_fp ); # endif /*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> * these routines written by Kevin Volk */ /* tnu array will be used to hold stellar array * * * NATLAS is the number of accessable models in the new model set * record 3 is the wavelength grid, record 2 holds the log(g) values, * and record 1 holds the T_eff values, so there are NATLAS+3 records * in the main file * * the values are all in integer format, usually of log10 quantities * real flux(ncell) */ /* these arrays hold the teff and log(g) values for the models in the * set; if the parameters are in range, interpolation is done first in * temperature and then in log(g) * * start by reading in the parameters * * * add together path to form DoAtlas file * * 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, "atlas.mod" ); } else { /* path not set, look here */ strcpy( chLine , "atlas.mod" ); } /* put junk in indices so we know if set */ for( j=0; j < 11; j++ ) { for( k=0; k < 61; k++ ) { jval[k][j] = -1; } } /* gjf open(unit=92,file='atlas.mod', */ if( lgTrace ) { fprintf( ioQQQ, " About to open file=%80.80s\n", chLine); } if( (ioIN = fopen( chLine , "rb" )) == NULL ) { fprintf( ioQQQ, "Error: Kurucz stellar atmosphere file atlas.mod not found.\n" ); fprintf( ioQQQ, " The path I tried:%80.80s\n", chLine ); fprintf( ioQQQ, " And here comes its octal representation:\n" ); fprintf( ioQQQ, " " ); for( i=0; i < 80; i++ ) { fprintf( ioQQQ, "%c%4x", chLine[i], chLine[i] ); } fprintf( ioQQQ, "\n" ); /* write(qq,'(1X,A1,O4)') (chLine(i:i),chLine(i:i),i=1,80) */ puts( "[Stop in DoAtlas]" ); exit(1); } if( lgTrace ) { fprintf( ioQQQ, " About to get wavelengths\n" ); } /* read in the saved cloudy energy scale so we can confirm this is a good image */ if( abs( fread( rfield.flux, 1, sizeof(rfield.flux),ioIN ) - sizeof(rfield.flux) ) !=0 ) { fprintf( ioQQQ, " problem trying to read atlas.mod wavelengths \n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(rfield.flux) ); puts( "[Stop in DoAtlas]" ); exit(1); } for( j=0; j < NCELL; j++ ) { /* tnu holds the frequency arrays */ IncidSpec.tnu[j][IncidSpec.nspec-1] = rfield.flux[j]; } if( lgTrace ) { fprintf( ioQQQ, " About to get temps\n" ); } if( abs( fread( rfield.flux, 1, sizeof(rfield.flux),ioIN ) - sizeof(rfield.flux) ) !=0 ) { fprintf( ioQQQ, " problem trying to read atlas.mod temperatures \n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(rfield.flux) ); puts( "[Stop in DoAtlas]" ); exit(1); } /* now copy over to actual teff array */ for( j=0; j < NATLAS; j++ ) { teff[j] = rfield.flux[j]; } if( lgTrace ) { fprintf( ioQQQ, " About to get gravs\n" ); } if( abs( fread( rfield.flux, 1, sizeof(rfield.flux),ioIN ) - sizeof(rfield.flux) ) !=0 ) { fprintf( ioQQQ, " problem trying to read atlas.mod gravities \n" ); fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n", (long)sizeof(rfield.flux) ); puts( "[Stop in DoAtlas]" ); exit(1); } for( j=0; j < NATLAS; j++ ) { xlogg[j] = rfield.flux[j]; } for( j=0; j < NATLAS; j++ ) { jj = (long)(2.*xlogg[j]+0.001) + 1; if( jj < 1 || jj > 11 ) { fprintf( ioQQQ, " DoAtlas finds insane jj:%5ld j, xlogg(j),=%5ld%10.2e\n", jj, j, xlogg[j] ); puts( "[Stop in DoAtlas]" ); exit(1); } for( k=0; k < 61; k++ ) { if( teff[j] == tval[k] ) jval[k][jj-1] = j + 1; } } /* determine which model or models are needed.*/ /* TEMP and ALOGG are desired temp and log g */ jmod = -1; for( j=0; j < NATLAS; ++j ) { dev1 = fabs(temp/teff[j]-1.); if( xlogg[j] > 0.4 ) { dev2 = fabs(alogg/xlogg[j]-1.); } else { dev2 = fabs(alogg-xlogg[j]); } /* this is because log(g) can be zero for the new set of models * this is test for exact model, result is jmod still -1 if hit */ if( dev1 <= 0.001 && dev2 <= 0.001 ) { /* +3 below is to skip first three records, which are cloudy energy * scale, teff, and log g */ jmod = j; break; } } if( alogg < 0. || alogg > 5. ) { /* this rejects temps outside the range of models */ fprintf( ioQQQ, " Gravity outside range of 0.0 to 5.0\n" ); puts( "[Stop in DoAtlas]" ); exit(1); } /* following blocks either a) interpolate or b) get exact model */ if( jmod == -1 ) { if( temp < 3500. || temp > 50000. ) { /* this rejects values outside the range of model temperatures */ fprintf( ioQQQ, " Temperature outside range of 3500 to 50,000K\n" ); puts( "[Stop in DoAtlas]" ); exit(1); } if( lgTrace ) { fprintf( ioQQQ, " About to call GETATL.\n" ); } /* DoAtlasInterp carries out the interpolation in t_eff and log(g) */ DoAtlasInterp(temp,alogg,jval,tval,&ierr,teff,xlogg , ioIN ); if( ierr == 1 ) { fprintf( ioQQQ, " GETATL returned ierr=1, stop.\n" ); puts( "[Stop in DoAtlas]" ); exit(1); } } else { if( lgTrace ) { fprintf( ioQQQ, " Exact model about to be read in.\n" ); } /* this branch, we hit a proper model, use exactly that one */ if( fseek(ioIN, (jmod+4)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error seeking exact Atlas atmosphere%4ld\n", jmod ); puts( "[Stop in DoAtlas]" ); exit(1); } if( abs( fread( rfield.flux, 1 , sizeof(rfield.flux), ioIN ) - sizeof(rfield.flux)) != 0 ) { fprintf( ioQQQ, " Error trying to read exact Atlas atmosphere%4ld\n", jmod ); puts( "[Stop in DoAtlas]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< K. Volk s Kurucz model %3ld" "read in. T_eff = %8.2f log(g) = %8.5f >>> *\n", jmod+1, teff[jmod], xlogg[jmod] ); } } for( j=0; j < NCELL; j++ ) { /* TSLOP is f_nu for each TNU point */ IncidSpec.tslop[j][IncidSpec.nspec-1] = rfield.flux[j]; } *nstar = NCELL; fclose(ioIN); # ifdef DEBUG_FUN fputs( " <->DoAtlas()\n", debug_fp ); # endif return; } /*DoAtlasInterp get one of the Atlas model atmospheres, coded by K. Volk */ void DoAtlasInterp(double temp, double alogg, long int jval[][11], float tval[], long int *ierr, float teff[], float xlogg[], FILE *ioIN ) { long int i, ipGrav, j, ipTeff, ipGravUp , ipTeffUp ; /* following used in binary read with fread, do not change to * double unless binary file changed too */ float flux1[NCELL], flux2[NCELL]; double fr1, fr2, fr3, fr4, xval; /* this is offset in call to fseek to skip initial records */ const long iOff = 3; # ifdef DEBUG_FUN fputs( "<+>DoAtlasInterp()\n", debug_fp ); # endif /* parameter (ndata=865,nmods=412) */ for( j=0; j < 61; j++ ) { for( i=9; i >= 0; i-- ) { if( jval[j][i] == -1 ) jval[j][i] = jval[j][i+1]; } } /* determine which models are to be interpolated */ ipTeff = -1; for( j=0; j < 60; j++ ) { if( tval[j] <= temp && tval[j+1] > temp ) { ipTeff = j; break; } } /* nb - ipTeff is pointer within c array */ if( ipTeff < 0 ) { if( tval[60] == temp ) { ipTeff = 60; } else { fprintf( ioQQQ, " Requested temperature of%10.2e is not within the range%10.2e%10.2e\n", temp, tval[0], tval[60] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } } ipTeffUp = MIN2( ipTeff+1 , 60 ); xval = 2.*alogg ; ipGrav = (long)(xval); ipGrav = MIN2( ipGrav , 10 ); ipGravUp = MIN2( ipGrav+1 , 10 ); /* Interpolation is carried out first in Teff and then in LOG(g). * Two Kurucz models are read in for the lower LOG(g) value first, which * I denote models A and B with B the higher T model [ipTeff+1 rather than ipTeff]. * unless alogg is 5.0 a second pair of Kurucz models are read in for the * higher LOG(g) value, denoted models C and D. Then the interpolation is * * RESULT = fr2*(fr3*B+fr4*A)+fr1*(fr3*D+fr4*C) * * when one has 4 models and * * RESULT = (fr3*B+fr4*A) * * when one has just 2 models. */ fr1 = xval - (float)(ipGrav); fr2 = 1. - fr1; if( ipTeffUp==ipTeff ) { fr3 = 0.; } else { fr3 = (temp - tval[ipTeff])/(tval[ipTeffUp] - tval[ipTeff]); } fr4 = 1. - fr3; /* this is check that we are within rational bounds */ if( fr3 < 0. || fr3 > 1. ) { fprintf( ioQQQ, "fr3 insanity in DoAtlasInterp\n" ); ShowMe(); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( fr2 < 0. || fr2 > 1. ) { fprintf( ioQQQ, "fr2 insanity in DoAtlasInterp\n" ); ShowMe(); puts( "[Stop in DoAtlasInterp]" ); exit(1); } /* at this point we should have valid gravity (ipGrav) and * teff (ipTeff) pointers */ if( fseek(ioIN, (jval[ipTeff][ipGrav]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error1 seeking Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( abs( fread( flux1, 1 , sizeof(flux1), ioIN ) - sizeof(flux1)) != 0 ) { fprintf( ioQQQ, " Error1 trying to read Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< K. Volk s Kurucz model %3ld" "read in. T_eff = %8.2f LOG(g) = %8.5f >>> *\n", jval[ipTeff][ipGrav], teff[jval[ipTeff][ipGrav]-2], xlogg[jval[ipTeff][ipGrav]-1] ); } if( fseek(ioIN, (jval[ipTeffUp][ipGrav]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error1 seeking Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( abs( fread( flux2, 1 , sizeof(flux2), ioIN ) - sizeof(flux2)) != 0 ) { fprintf( ioQQQ, " Error2 trying to read Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< K. Volk s Kurucz model %3ld" "read in. T_eff = %8.2f LOG(g) = %8.5f >>> *\n", jval[ipTeffUp][ipGrav], teff[jval[ipTeffUp][ipGrav]-2], xlogg[jval[ipTeffUp][ipGrav]-1] ); } for( j=0; j < NCELL; j++ ) { rfield.flux[j] = (float)( fr4*log10(MAX2(1e-37,flux1[j]) ) + fr3*log10(MAX2(1e-37,flux2[j]) ) ); } if( ipGrav < 10 ) { /* for LOG(g) < 5.0 repleat the interpolation at ht elarger LOG(g) then * carry out the interpolation in log(g) * >>chng 97 july 24, get rid of go to */ if( fseek(ioIN, (jval[ipTeff][ipGravUp]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error3 seeking Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( abs( fread( flux1, 1 , sizeof(flux1), ioIN ) - sizeof(flux1)) != 0 ) { fprintf( ioQQQ, " Error3 trying to read Atlas atmosphere%4ld\n", jval[ipTeff][ipGravUp] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< K. Volk s Kurucz model %3ld" "read in. T_eff = %8.2f log(g) = %8.5f >>> *\n", jval[ipTeff][ipGravUp], teff[jval[ipTeff][ipGravUp]-2], xlogg[jval[ipTeff][ipGravUp]-1] ); } if( fseek(ioIN, (jval[ipTeffUp][ipGravUp]+iOff)*sizeof(rfield.flux), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error4 seeking Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( abs( fread( flux2, 1 , sizeof(flux2), ioIN ) - sizeof(flux2)) != 0 ) { fprintf( ioQQQ, " Error4 trying to read Atlas atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoAtlasInterp]" ); exit(1); } if( called.lgTalk ) { fprintf( ioQQQ, " * <<< K. Volk s Kurucz model %3ld" "read in. T_eff = %8.2f log(g) = %8.5f >>> *\n", jval[ipTeffUp][ipGravUp], teff[jval[ipTeffUp][ipGravUp]-2], xlogg[jval[ipTeffUp][ipGravUp]-1] ); } for( j=0; j < NCELL; j++ ) { flux1[j] = (float)( fr4*log10(MAX2(1e-37,flux1[j] ) ) + fr3*log10(MAX2(1e-37,flux2[j] ) ) ); } for( j=0; j < NCELL; j++ ) { /* >>chng 96 dec 27, as per K Volk change, * to interpolation in log of space*/ if( rfield.flux[j]>-36. && flux1[j]>-36. ) { rfield.flux[j] = (float)pow(10.,fr2*rfield.flux[j] + fr1*flux1[j]); } else { rfield.flux[j] = 0.; } } } else { /* One gets here if log(g) is 5.0 and there is no higher log(g) pair of * models to interpolate. * */ for( j=0; j < NCELL; j++ ) { rfield.flux[j] = (float)pow(10.,rfield.flux[j]); if( rfield.flux[j]<1e-36 ) rfield.flux[j] = 0.; } } *ierr = 0; # ifdef DEBUG_FUN fputs( " <->DoAtlasInterp()\n", debug_fp ); # endif return; }