/*ConcatRauch K Volk's routine to compact the 66 Rauch atmospheres into a large ascii file.*/ /*RebinRauch rebin the rauch continuum grid onto the cloudy grid by K Volk */ /*DoRauch get one of the Rauch model atmospheres, coded by K. Volk */ /*GetRauch rebin Rauch stellar models to match energy grid of code*/ #include #include "cddefines.h" #include "opac.h" #include "printe82.h" #include "incidspec.h" #include "path.h" #include "called.h" #include "rfield.h" #include "showme.h" #include "rauch.h" /* this is the number of points in each Rauch stellar continuum */ #define NRAUCH 19951 /* number of atmospheres */ #ifdef NMODS #undef NMODS #endif #define NMODS 66 /* rebin the continuum, called by GetRauch below */ void RebinRauch( /* 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, float *slope); /*========================================================================*/ void ConcatRauch(void) { char chFName[21], chLine[133], /* used for getting input lines */ chFileName[133] , /* file name for star input atmos */ chSuff[11]; /* the suffix for the star file */ /* As the wavelengths are always the same only the fluxes * will be written to the output file. * * Before running this program issue the following command where the Rauch * model atmosphere files are kept (050000_50.wf and so on) * * ls *.wf > rauchmods.list * * and check to see that there are 66 lines in the file. The file name * "1000000_90_solar.wf" then has to be moved to the end of the file for * the models to be written to the ascii output file in the proper order. */ /* concatenated output will go to this file */ FILE *ioOut, /* pointer to output file we came here to create*/ *ioIn ; /* pointer to input files we will reaed */ long int i, j, k; float wl; static char chFile[NMODS][11]= {"050000_50","050000_60","050000_70", "050000_80","060000_50","060000_60","060000_70","060000_80", "070000_50","070000_60","070000_70","070000_80","080000_50", "080000_60","080000_70","080000_80","090000_50","090000_60", "090000_70","090000_80","100000_50","100000_60","100000_70", "100000_80","110000_60","110000_70","110000_80","120000_60", "120000_70","120000_80","130000_60","130000_70","130000_80", "140000_60","140000_70","140000_80","150000_60","150000_70", "150000_80","160000_60","160000_70","160000_80","170000_60", "170000_70","170000_80","180000_60","180000_70","180000_80", "190000_60","190000_70","190000_80","200000_70","200000_80", "200000_90","300000_70","300000_80","300000_90","400000_80", "400000_90","500000_80","500000_90","600000_90","700000_90", "800000_90","900000_90","1000000_90"}; # ifdef DEBUG_FUN fputs( "<+>ConcatRauch()\n", debug_fp ); # endif fprintf( ioQQQ, "ConcatRauch on the job. This will take a while.\n" ); /* loop over two types of stellar atmospheres, solar and halo */ for( k=0; k < 2; k++ ) { if( k == 0 ) { strcpy( chFName, "rauch_solar.ascii" ); strcpy( chSuff, "_solar.wf" ); fprintf( ioQQQ, "Doing solar abundance set....\n" ); } else { strcpy( chFName, "rauch_halo.ascii" ); strcpy( chSuff, "_halo.wf" ); fprintf( ioQQQ, "Doing halo abundance set....\n" ); } ioOut = fopen( chFName , "w" ); /* check that we got it ok */ if( ioOut == NULL ) { fprintf(ioQQQ,"error opening output file for Rauch atmospheres.\n"); fprintf(ioQQQ,"This is routine concatrauch.\n"); puts( "[Stop in concatrauch]" ); exit(1); } for( i=0; i < NMODS; i++ ) { /* must create name of next stellar atmosphere */ strcpy( chFileName , chFile[i] ); strcat( chFileName , chSuff ); /* now open next stellar atmosphere for reading*/ ioIn = fopen( chFileName , "r" ); if( ioIn == NULL ) { fprintf(ioQQQ,"error opening input file for Rauch atmospheres.\n"); fprintf(ioQQQ,"name was %s\n",chFileName); fprintf(ioQQQ,"This is routine concatrauch.\n"); puts( "[Stop in concatrauch]" ); exit(1); } /* print to screen to show progress */ fprintf( ioQQQ, " star%3li %s\n", i, chFile[i]); /* throw away the header information */ for( j=0; j < 33; j++ ) { if(fgets( chLine , sizeof(chLine) , ioIn )==NULL ) { fprintf( ioQQQ, " ConcatRauch error in header atmosphile file %4ld%4ld\n", i, j ); puts( "[Stop in concatrauch]" ); exit(1); } } for( j=0; j < 19951; j++ ) { /* get the input line */ if(fgets( chLine , sizeof(chLine) , ioIn )==NULL ) { fprintf( ioQQQ, " ConcatRauch error in header atmosphile file %4ld%4ld\n", i, j ); puts( "[Stop in concatrauch]" ); exit(1); } /* scan off wavelength (don't save) and flux)*/ sscanf( chLine , "%6f%10e", &wl, &opac.OpacStack[j] ); } /* finished - close the unit */ fclose(ioIn); /* now write to output file */ PrintE93( ioOut, opac.OpacStack[0] ); for( j=1; j < 19951; j++ ) { PrintE93( ioOut, opac.OpacStack[j] ); /* want to have 9 numbers per line */ if(!((j+1)%9)) fprintf(ioOut,"\n"); } /* need to throw a newline if we did not end on an exact line */ if(((j+1)%9)) fprintf(ioOut,"\n"); } fclose(ioOut); } puts( "[Stop in concatrauch]" ); exit(1); } /*GetRauch rebin Rauch stellar models to match energy grid of code * limit to number of cells in frequency array * NB is this every is changed, also need to change formats in * Kevin Volk's stellar atmospheres codes */ /* helper routine located below */ /*========================================================================*/ void GetRauch(void) { FILE /* used for input */ * ioIN , /* used for output */ * ioOUT ; /* these will be malloced into large work arrays*/ float *StarEner , *StarFlux , *CloudyFlux , *scratch ; /* used to save line image */ char chLine[100]; long int i, imod, j, k; /* these will be malloced into large work arrays*/ # ifdef DEBUG_FUN fputs( "<+>GetRauch()\n", debug_fp ); # endif /* malloc some workspace */ StarEner = (float *)malloc( sizeof(float)*NRAUCH ); if( StarEner == NULL ) { printf( " not enough memory to allocate StarEner in GetRauch\n" ); puts( "[Stop in GetRauch]" ); exit(1); } /* following has extra 9 since we must read last line, which can be partly empty*/ scratch = (float *)malloc( sizeof(float)*(NRAUCH+9) ); if( scratch == NULL ) { printf( " not enough memory to allocate scratch in GetRauch\n" ); puts( "[Stop in GetRauch]" ); exit(1); } StarFlux = (float *)malloc( sizeof(float)*NRAUCH ); if( StarFlux == NULL ) { printf( " not enough memory to allocate StarFlux in GetRauch\n" ); puts( "[Stop in GetRauch]" ); exit(1); } CloudyFlux = (float *)malloc( sizeof(float)*NCELL ); if( CloudyFlux == NULL ) { printf( " not enough memory to allocate CloudyFlux in GetRauch\n" ); puts( "[Stop in GetRauch]" ); exit(1); } /* This subroutine re-bins the Rauch atmosphere fluxes to match the Cloudy * grid. The wavelengths of the Rauch models are 5.0 to 2000.0 angstroms * with a step size of 0.1 angstroms, 19951 points each. Thus the wavelength * grid is not saved in the file of ascii values. The values in the file are * F_lambda values in erg/s/cm^2/cm. They have to be multiplied by pi to get * the "real" surface flux. Then I convert them into F_nu in erg/s/cm^2/Hz * by taking 1.e-18 * wavelength^2 * F_lambda / c. The factor 1.e-18 is * due to the wavelengths being in angstroms.*/ /* With some many points I need to average over the Cloudy grid rather than * interpolating as is the case for the Kurucz or Werner models. At short * wavelengths the two grids are of similar spacing. */ /* nvals is the number of values per Rauch model. nmods is the total * number of models in the set.*/ fprintf( ioQQQ, " GetRauch on the job.\n" ); fprintf( ioQQQ, " There are a total of %i models, and I will need to do this twice.\n",NMODS ); /* StarEner is the Rauch stellar atmosphere wavelength grid*/ for( i=0; i < NRAUCH; i++ ) { /* notice that this is in decreasing wavelength order */ StarFlux[i] = (float)(2000. - 0.1*(float)(i)); } /* convert continuum over to increasing energy in Ryd */ for( i=0; i < NRAUCH; i++ ) { StarEner[i] = (float)(1e10/(StarFlux[i]*10967758.30)); } for( k=1; k <= 2; k++ ) { if( k == 1 ) { /* this is the solar rauch grid */ /*if( (ioIN = fopen( "rauch_solar.ascii", "r" ) ) == NULL )*/ if( (ioIN = fopen( "rauch_solar.ascii", "r" ) ) == NULL ) { fprintf( ioQQQ, " GetRauch could not find rauch_solar.ascii.\n" ); puts( "[Stop in getrauch]" ); exit(1); } fprintf( ioQQQ, " GetRauch got rauch_solar.ascii.\n" ); /* the frequencies are in increasing order; when the fluxes are * read in they will have to be put in reverse order since they go from * 5 to 200 angstroms. */ if( (ioOUT = fopen( "rauch_solar.mod", "wb" ) ) == NULL ) { fprintf( ioQQQ, " GetRauch failed creating rauch_solar.mod.\n" ); puts( "[Stop in getrauch]" ); exit(1); } } else { /* this is the halo rauch grid */ if( (ioIN = fopen( "rauch_halo.ascii", "r" ) ) == NULL ) { fprintf( ioQQQ, " GetRauch could not find rauch_halo.ascii.\n" ); puts( "[Stop in getrauch]" ); exit(1); } fprintf( ioQQQ, " GetRauch got rauch_halo.ascii.\n" ); if( (ioOUT = fopen( "rauch_halo.mod", "wb" ) ) == NULL ) { fprintf( ioQQQ, " GetRauch failed creating rauch_halo.mod.\n" ); puts( "[Stop in getrauch]" ); exit(1); } } /* The Rauch models go down to 5 Angstroms (about 182 Rydb) but I choose to * terminate the array at 150 Rydb since the fluxes are extremely small at * the shortest wavelengths -- 17 orders of magnitude below peak for the * 500000 K modelT. This differs a bit from the handling of the Werner models * even though the temperature range is similar.*/ /* 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, " GetRauch failed writng anu array.\n" ); puts( "[Stop in getrauch]" ); 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). */ for( imod=0; imod < NMODS; imod++ ) { /* now read the energy grid */ i = 0; while( i < NRAUCH ) { if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL ) { fprintf( ioQQQ, " GetRauch failed reading star flux.\n" ); puts( "[Stop in getrauch]" ); exit(1); } sscanf( chLine , "%9f %9f %9f %9f %9f %9f %9f %9f %9f" , &scratch[i ], &scratch[i+1], &scratch[i+2], &scratch[i+3], &scratch[i+4], &scratch[i+5], &scratch[i+6], &scratch[i+7], &scratch[i+8]); /* increment counter, nine numbers per line*/ i += 9; } for( j=0; j < NRAUCH; j++ ) { /* flip values around (wavelength to frequency) and convert to f_nu */ StarFlux[j] = scratch[NRAUCH-1-j]/(StarEner[j]*StarEner[j]); } /* the re-binned values are returned in the "CloudyFlux" array */ RebinRauch(StarEner,StarFlux,CloudyFlux,scratch); /*{ FILE *ioBUG; ioBUG = fopen("test.txt","w"); for( i=0; iGetRauch()\n", debug_fp ); # endif return; } /*RebinRauch rebin the rauch continuum grid onto the cloudy grid by K Volk */ /*========================================================================*/ void RebinRauch( /* 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, float *slope) { long int i, ipHi, ipLo, j; double EdgeHigh, EdgeLow, slope2, sum, v1, v2, value, x1, x2; # ifdef DEBUG_FUN fputs( "<+>RebinRauch()\n", debug_fp ); # endif /* Arrays passed in from "GetRauch" */ /* This is mostly a copy of "RebinAtlas" except that I pass in the arrays * rather than the integer values, since both values needed will * be fixed before this routine is called and since the model arrays are * much larger than the Cloudy grid size. Also I do linear interpolation * rather than log-log interpolation here. * * StarFlux holds the current model flux densities F_nu in erg/s/cm^2/Hz. * wlin holds the wavelengths in Angstroms. * CloudyFlux will have the rebinned fluxes in it upon return. * * limit to number of cells in frequency array * NB is this every is changed, also need to change formats in * Kevin Volk's stellar atmospheres codes */ /* 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. * * Convert wavelengths from angstroms to rydberg; in the calling routine * I assured that the wavelengths are decreasing in the array.*/ /* Calculate the slopes for interpolation between the Rauch model points. * Since the model points are closely spaced I will use linear interpolation.*/ for( j=0; j < NRAUCH-1; j++ ) { slope[j] = (StarFlux[j+1] - StarFlux[j])/(StarEner[j+1] - StarEner[j]); } /* this is the Rayleigh-Jeans slope in F_nu, *for "low" frequency extrapolation.*/ slope2 = 2.0; /* Actually the extrapolation does not match too well for the lower temperature * models, but if the long wavelength stellar continuum level matters in the * results then something strange is going on anyway.*/ for( j=0; j StarEner[NRAUCH-1] )*/ if( EdgeHigh > StarEner[NRAUCH-1] ) { /* cloudy continuum above stellar continuum, set to zero */ CloudyFlux[j] = 0.0; } else { /* find lower energy bound */ ipLo = 0; for( i=0; i < (NRAUCH - 1); i++ ) { /* find where stellar continuum maps onto cloudy continuum */ if( StarEner[i] <= EdgeLow && StarEner[i+1] > EdgeLow ) { ipLo = i; break; } } /* now find upper energy bound */ ipHi = ipLo; for( i=MAX2(0,ipLo-1); i < (NRAUCH - 1); i++ ) { if( StarEner[i] <= EdgeHigh && StarEner[i+1] > EdgeHigh ) { ipHi = i; break; } } if( ipHi == ipLo || ipLo>ipHi) { /* case where cloudy cell bigger than stellar cell */ CloudyFlux[j] = StarFlux[ipHi-1] + slope[ipHi-1]*(rfield.anu[j] - StarEner[ipHi-1]); } else { sum = 0.; for( i=ipLo; i <= ipHi; i++ ) { v1 = MAX2(StarEner[i],EdgeLow); v2 = MIN2(StarEner[i+1],EdgeHigh); value = v2 - v1; x1 = StarFlux[i] + slope[i]*(v1 - StarEner[i]); x2 = StarFlux[i] + slope[i]*(v2 - StarEner[i]); sum += (x1 + x2)*value; } sum /= 2.*rfield.widflx[j]; CloudyFlux[j] = (float)sum; } } } } # ifdef DEBUG_FUN fputs( " <->RebinRauch()\n", debug_fp ); # endif return; } /*DoRauch get one of the Rauch model atmospheres, coded by K. Volk */ /*lint -e713 unsigned to signed */ /*lint -e737 signed to unsigned */ /*========================================================================*/ void DoRauch(long int *nstar, double temp, double alogg, long int iabund) { /* pointer to the rebinned binary grids */ FILE *ioIN ; char chLine[81]; long int i, ipGrav, ipGravUp, j, ipTeff , ipTeffUp; /* these are used to get the binary files, * they must not be changed to double without change those files too */ float flux1[NCELL], flux2[NCELL]; double fact1, fact2, fr1, fr2, fr3, fr4, xval; static double teff[NMODS]={50000.,50000.,50000.,50000.,60000.,60000., 60000.,60000.,70000.,70000.,70000.,70000.,80000.,80000.,80000., 80000.,90000.,90000.,90000.,90000.,100000.,100000.,100000.,100000., 110000.,110000.,110000.,120000.,120000.,120000.,130000.,130000., 130000.,140000.,140000.,140000.,150000.,150000.,150000.,160000., 160000.,160000.,170000.,170000.,170000.,180000.,180000.,180000., 190000.,190000.,190000.,200000.,200000.,200000.,300000.,300000., 300000.,400000.,400000.,500000.,500000.,600000.,700000.,800000., 900000.,1000000.}; static double xlogg[NMODS]={5.,6.,7.,8.,5.,6.,7.,8.,5.,6.,7.,8., 5.,6.,7.,8.,5.,6.,7.,8.,5.,6.,7.,8.,6.,7.,8.,6.,7.,8.,6.,7., 8.,6.,7.,8.,6.,7.,8.,6.,7.,8.,6.,7.,8.,6.,7.,8.,6.,7.,8.,7., 8.,9.,7.,8.,9.,8.,9.,8.,9.,9.,9.,9.,9.,9.}; static double tval[24]={50000.,60000.,70000.,80000.,90000.,100000., 110000.,120000.,130000.,140000.,150000.,160000.,170000.,180000., 190000.,200000.,300000.,400000.,500000.,600000.,700000.,800000., 900000.,1000000.}; /* The 0th offset block is the wavelength scale, * the first model is number 2 in jval, so you must subtract * 1 when using jval to point to tval or xlogg */ static long jval[24][5]={ {2,3,4,5,5}, {6,7,8,9,9}, {10,11,12,13,13}, {14,15,16,17,17}, {18,19,20,21,21}, {22,23,24,25,25}, {26,26,27,28,28}, {29,29,30,31,31}, {32,32,33,34,34}, {35,35,36,37,37}, {38,38,39,40,40}, {41,41,42,43,43}, {44,44,45,46,46}, {47,47,48,49,49}, {50,50,51,52,52}, {53,53,53,54,55}, {56,56,56,57,58}, {59,59,59,59,60}, {61,61,61,61,62}, {63,63,63,63,63}, {64,64,64,64,64}, {65,65,65,65,65}, {66,66,66,66,66}, {67,67,67,67,67} }; # ifdef DEBUG_FUN fputs( "<+>DoRauch()\n", debug_fp ); # endif /* get one of the Rauch model atmospheres, coded by K. Volk */ /* limit to number of cells in frequency array * NB is this every is changed, also need to change formats in * Kevin Volk's stellar atmospheres codes */ if( lgDataPathSet ) { /* path is set, generate full path name with file */ strcpy( chLine , chDataPath ); if( iabund == 0 ) { /* solar abundances, the default */ strcat( chLine, "rauch_solar.mod" ); } else { /* halo abundances, set with halo keyword */ strcat( chLine, "rauch_halo.mod" ); } } else { /* path not set, look here */ if( iabund == 0 ) { strcpy( chLine, "rauch_solar.mod" ); } else { strcpy( chLine, "rauch_halo.mod" ); } } if( (ioIN = fopen( chLine , "rb" )) == NULL ) { /* something went wrong */ fprintf( ioQQQ, "Error: Rauch stellar atmosphere file rauch_solar.mod or rauch_halo.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" ); puts( "[Stop in rauchmod]" ); exit(1); } /* read in the saved cloudy energy scale so we can confirm this is a good image */ if( abs( fread( flux1, 1, sizeof(flux1),ioIN ) - sizeof(flux1) ) !=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(flux1) ); puts( "[Stop in DoRauch]" ); exit(1); } for( j=0; j < NCELL; j++ ) { /* save energy scale for check against code's in conorm * (scale not defined when this routine called */ IncidSpec.tnu[j][IncidSpec.nspec-1] = (float)flux1[j]; } /* determine which models are to be interpolated */ ipTeff = -1; for( j=0; j < 23; j++ ) { if( tval[j] <= temp && tval[j+1] > temp ) { /* this points to model at or below desired temperature */ ipTeff = j; break; } } /* this happens if above loop never got a hit */ if( ipTeff<0 ) { /* above loop could not catch highest temperature */ if( tval[23] == temp ) { ipTeff = 23; } else { fprintf( ioQQQ, " Requested temperature of%11.2e is not within the range" "%10.2e to e%10.2e\n", temp, tval[0], tval[23] ); puts( "[Stop in DoRauch]" ); exit(1); } } /* we now have fully validated pointers to lower and upper limits to * stellar temperature scale */ ipTeffUp = MIN2( ipTeff+1 , 23 ); /* now generate same types of pointers within gravity scale, * which has 1 dex steps, log 5 = 5, 6, 7, and 8 for cooler stars, * hottest stars only 8 */ xval = (alogg - 4.0); /* >>>chng 99 apr 28, from above to below, apparently left on Fort scale, * move onto c scale for addressing array */ xval = (alogg - 5.0); ipGrav = (long)(xval); ipGrav = MAX2(0, ipGrav); ipGrav = MIN2( ipGrav , 4 ); ipGravUp = MIN2( ipGrav+1, 4); /* Interpolation is carried out first in Teff and then in log(g). * Two Rauch 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 8.0 a second pair of Rauch 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. */ /* this sets up interpolation fractions in log g */ fr1 = xval - (float)(ipGrav); fr2 = (1. - fr1); /* now set up interplation fractions in Teff */ if( ipTeffUp == ipTeff ) { fr3 = 0.; } else { fr3 = ((temp - tval[ipTeff])/(tval[ipTeffUp] - tval[ipTeff])); } fr4 = 1. - fr3; if( fr3>1. || fr3 < 0. ) { fprintf( ioQQQ, " fr3 insanity in DoRauch\n" ); ShowMe(); puts( "[Stop in DoRauch]" ); exit(1); } if( fr2>1. || fr2 < 0. ) { fprintf( ioQQQ, " fr2 insanity in DoRauch\n" ); ShowMe(); puts( "[Stop in DoRauch]" ); exit(1); } /* skip over jval-1 stars */ if( fseek(ioIN, (jval[ipTeff][ipGrav]-1)*sizeof(flux1), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error1 seeking Rauch atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } if( abs( fread( flux1, 1 , sizeof(flux1), ioIN ) - sizeof(flux1)) != 0 ) { fprintf( ioQQQ, " Error1 trying to read Rauch atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } if( called.lgTalk ) { /* say we read in the first atmosphere */ fprintf( ioQQQ, " *<<< Rauch model %3ld READ IN" ". T_eff = %10.2f log(g) = %8.5f >>> *\n", /* -1 below to be same as C90 */ jval[ipTeff][ipGrav]-1, teff[jval[ipTeff][ipGrav]-2], xlogg[jval[ipTeff][ipGrav]-2] ); } if( fseek(ioIN, (jval[ipTeffUp][ipGrav]-1)*sizeof(flux2), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error2 seeking Rauch atmosphere%4ld\n", jval[ipTeffUp][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } if( abs( fread( flux2, 1 , sizeof(flux2), ioIN ) - sizeof(flux2)) != 0 ) { fprintf( ioQQQ, " Error2 trying to read Rauch atmosphere%4ld\n", jval[ipTeffUp][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } /* convert to logs since we will interpolate in log flux */ for( i=0; i>> *\n", jval[ipTeffUp][ipGrav]-1, teff[jval[ipTeffUp][ipGrav]-2], xlogg[jval[ipTeffUp][ipGrav]-2] ); } /* The following is an approximate scaling to use for the range of * temperatures above 200000 K where the temperature steps are large * and thus the interpolations are over large ranges. For the lower * temperatures I assume that there is no need for this since the * step of 10000 K is not so large. * * It should be remembered that this interpolation is not exact, and * the possible error at high temperatures might be large enough to * matter. */ if( temp > 200000. ) { fact1 = log10(temp/tval[ipTeff])*4.; fact2 = log10(temp/tval[ipTeffUp])*4.; } else { fact1 = 0.; fact2 = 0.; } /* save interpolated continuum here for now, will use again below */ for( j=0; j < NCELL; j++ ) { rfield.flux[j] = (float)(fr4*(flux1[j] + fact1) + fr3*(flux2[j] + fact2)); } /* ipGrav is offset from lowest gravity, so ranges from 0 through 3 */ if( ipGrav < 4 ) { /* for log(g) < 8.0 repeat the interpolation at the larger log(g) then * carry out the interpolation in log(g) */ if( fseek(ioIN, (jval[ipTeff][ipGravUp]-1)*sizeof(flux1), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error3 seeking Rauch atmosphere%4ld\n", jval[ipTeff][ipGravUp] ); puts( "[Stop in DoRauch]" ); exit(1); } if( abs( fread( flux1, 1 , sizeof(flux1), ioIN ) - sizeof(flux1)) != 0 ) { fprintf( ioQQQ, " Error3 trying to read Rauch atmosphere%4ld\n", jval[ipTeff][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } /* convert to logs since we will interpolate in log flux */ for( i=0; i>> *\n", jval[ipTeff][ipGravUp]-1, teff[jval[ipTeff][ipGravUp]-2], xlogg[jval[ipTeff][ipGravUp]-2] ); } if( fseek(ioIN, (jval[ipTeffUp][ipGravUp]-1)*sizeof(flux2), SEEK_SET )!= 0 ) { fprintf( ioQQQ, " Error2 seeking Rauch atmosphere%4ld\n", jval[ipTeffUp][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } if( abs( fread( flux2, 1 , sizeof(flux2), ioIN ) - sizeof(flux2)) != 0 ) { fprintf( ioQQQ, " Error2 trying to read Rauch atmosphere%4ld\n", jval[ipTeffUp][ipGrav] ); puts( "[Stop in DoRauch]" ); exit(1); } /* convert to logs since we will interpolate in log flux */ for( i=0; i>> *\n", jval[ipTeffUp][ipGravUp]-1, teff[jval[ipTeffUp][ipGravUp]-2], xlogg[jval[ipTeffUp][ipGravUp]-2] ); } /* at this stage things are still log flux, so interpolate */ for( j=0; j < NCELL; j++ ) { flux1[j] = (float)((flux1[j] + fact1)*fr4 + (flux2[j] + fact2)*fr3); } /* now convert back to linear, set to zero if very small */ for( j=0; j < NCELL; j++ ) { rfield.flux[j] = (float)pow(10.,fr2*rfield.flux[j] + fr1*flux1[j]); if( rfield.flux[j] < 1e-37 ) rfield.flux[j] = 0.; } } else { /* One gets here if log(g) is 8.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-37 ) rfield.flux[j] = 0.; } } /* Note on the interpolation: 15 August 1997 (Kevin Volk) * * I tested the total flux from the interpolated grids. I found that the * ratio of actual to expected total flux is about 1.05 for 55000 K and * log(g) of 7.5, between 0.99 and 1.01 for temperature values between * 60000 K and 200000 K, and about 1.15 for 250000 K declining to 0.95 * for 500000 K. So for most of the range things are O.K. but where * the interpolation interval is large and at the highest temperatures * there is a noticeable scaling. * * I do not know why the hottest models require a scaling factor of less than * 1.0. This may mean that the extrapolation beyond 2000 Angstroms is bad * although one would think that at such high temperatures it would not * matter. * * Since Cloudy checks the scaling elsewhere there is no need to re-scale * things here, but this inaccuracy should be kept in mind. * */ 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 );; /*lint +e713 unsigned to signed */ /*lint +e737 signed to unsigned */ # ifdef DEBUG_FUN fputs( " <->DoRauch()\n", debug_fp ); # endif return; }