/******************************** * * FILE INPUT.C * * All the input functions are in this file. * ************************************/ #include "header.h" void readParams( void ) /******************** * * This function reads the input parameters from a file named "starpars". * The function assumes a FITS-like data file with keywords, then blanks, * then a set of parameters. Any extra parameters are ignored. * The active keywords are: * * Control keywords: * STAR1= ON or OFF * STAR1SPOTS= ON or OFF * STAR2= ON or OFF * STAR2SPOTS= ON or OFF * DISK= ON or OFF * THIRDLIGHT= ON or OFF * REFLECTION= ON or OFF * WRITEGRIDS= ON or OFF * a diagnostic tool. * RUNTYPE= LIGHTCURVE or * LINEIMAGE plotprogram * POINTIMAGE plotprogram * where "plotprogram" is the ascii name of the * ploting program, such as, "supermongo" * NORMALIZE= OFF * Do not normalize the light curve. * SETVALUE phase value * Normalize the flux in the synthetic light * curve by setting the flux to "value" at * phase "phase". * MAXVALUE maxvalue * Normalize the flux to maxvalue at the * phase of maximum flux. * LSQR phase1 phase2 * Normalize by an amount that minimizes the * variance of the fit to the observed light * curve. The fit will be made over phases * from "phase1" to "phase2". * PHASESOURCE= EVENSPACE firstphase lastphase deltaphase It is okay to have firstphase = lastphase in order to get a single phase point. * FILE phasefilename * See the comment at the beginning of function * OrbPhases() for the meanings of parameters. * READDATA= OFF or * ON datafilename * See the comment at the beginning of function * readData() for the file format. * COMMENT= Comment line. * blank line Blank lines have no effect. * END Terminates reading of the parameter file. * * Keywords referring to the whole system: * Q= q * Mass ratio m1 / m2 * MASS1= mass1 * Mass of star 1. Only needed for radial * velocity and line profile calculations. * INCLINATION= inclination * Orbital inclination in degrees. * FILTERWAVEL= filtername wavelength * The filter name should be something like: * F135W, U, B, V, R, I, J, H, K, etc. * The filter name is only used if the fluxes * come from the ITABLE option. * The wavelength is the effective wavelength * of the light curve. It is used if the fluxes * are calculated from the Planck function. * PHASEOFFSET= phaseoffset * The flux reported by the program for phase * phi will be calculated at phase * (phi - phaseoffset). * * Keywords referring to star 2 (with obvious equivalents for star 1). * THETAPOINTS2= thetapoints * Number of grid points in the theta direction * PHIPOINTS2= phipoints * Number of grid points in the phi direction * FRACRADIUS2= fracRadious * R2 / RL2; that is, the radius of star 2 * divided by the distance of the inner * Lagrangian point from the center of R2. * R2 is measured along the line joining the * centers of the stars. * POLETEMP2= poletemp * Temperature at the pole * GRAVBETA2= gravbeta * Gravity darkening beta: T ~ pow( g, beta) * FLUXSOURCE2= BB bblimbdark * Use Black Body fluxes calculated at the * effective wavelength of FILTER and use * bblimbdark for the linear limb darkening * coefficient where * limbdark = norm*(1 - u*(1-mu)) * ITABLE log(g) * Use the flux table itable.dat for the fluxes. * ALBEDO2= albedo * (actually 1-albedo) * STAR2SPOT= theta phi radius spotToverStarT * (angles in degrees; multiple spots allowed) * NOTE: theta = 0 points at the north pole * of star 2. * theta = 90, phi = 0 points at the * inner Lagrangian point * * Keywords referring to the ENTIRE the disk. * Note that the disk is always around star 1. * RPOINTSDISK= rpoints * Number of grid points in the r direction * Note: 3 grid points are on the edge of the * disk and the rest on the top of the disk. * PHIPOINTSDISK= phipoints * Number of disk grid points in the phi direction. * RDISKIN= innerR * RDISKOUT= outerR * The inner and outer radii of the disk. * The radii are in units of RL1, that is, the * the radii are innerR/RL1 and outerR/RL1. * Also, these are the radii of the disk * projected to the plane of the orbit, not * along the face of the disk. * If star 1 is ON, and if innerR is less * than the maximum radius of star 1, innerR * is set equal to the maximum radius of star 1. * FLAREANGLE= flareangle * Half-opening-angle of the disk in degrees. * DISKTRANSMIT= transmission * The transmission of the disk. * 1.0 = fully transparent. * 0.0 = opaque * * Keywords referring to the TOP of the disk. * TOPTEMP= POLYNOMIAL a0 a1 a2 a3 * Disk temperature is given by the polynomial * T = a0 + a1*r + a2*r^2 + a3*r^3 * where r is the radius projected to the * plane of the orbit. * STEADYSTATE To beta * Disk temperature is equal to the theoretical * distribution for an optically-thick, * steady-state disk. * FILE filename * The disk tempertures are read from a file * with name "filename". * TOPFLUXSOURCE= BB bblimbdark * Use Black Body fluxes calculated at the * effective wavelength of FILTER and use * bblimbdark for the linear limb darkening * coefficient where * limbdark = norm*(1 - u*(1-mu)) * ITABLE log(g) * Use the flux table itable.dat for the fluxes. * TOPSPOTTEMP= temperature * This puts a spot on top of the disk. * Turn off the spot by omitting this keyword * or by setting T < 0.0. * TOPSPOTR= rin rout * rin and rout are measured in the plane of * the orbit. If rin is less than the inner * radius of the disk it is set equal to the * inner radius. If rout is greater than the * outer radius of the disk it is set equal to * the outer radius. The spot is fit onto the * grid points by rounding its boundaries to * values at grid lines. * TOPSPOTPHI= phi1 phi2 * phi1 and phi2 are in degrees. * If phi2 is less than phi1, the spot goes * from phi1 through 360 degrees to phi2. * * Keywords referring to the RIM of the disk. * RIMTEMP= temperature * Turn off the spot by omitting this keyword * or by setting T < 0.0. * RIMFLUXSOURCE= BB bblimbdark * Use Black Body fluxes calculated at the * effective wavelength of FILTER and use * bblimbdark for the linear limb darkening * coefficient where * limbdark = norm*(1 - u*(1-mu)) * ITABLE log(g) * Use the flux table itable.dat for the fluxes. * RIMSPOTTEMP= temperature * This puts a spot on top of the disk. * Turn off the spot by omitting this keyword * or by setting T < 0.0. * RIMSPOTPHI= phi1 phi2 * If phi2 is less than phi1, the spot goes * from phi1 through 360 degrees to phi2. * The spot is fit to phi grid points by rounding * its boundaries to values at grid points. * * Keywords referring to the "third light" * FRACTION3= fraction * The fraction of the total flux that comes * from a constant background source. * PHASE3= phase * The orbital phase at which the fraction of the * flux coming from the third light is calculated. * *********************/ { char filename[40], inputline[81], keyword[20], param1[20], param2[20], param3[20], param4[20], param5[20]; double x, firstphase, lastphase, deltaphase; FILE *in; long int nfields; strcpy(filename, "starpars"); if( ( in = fopen(filename, "r")) == NULL ) quit("Cannot open file starpars."); while( fgets( inputline, 80, in) != NULL ) { /*********************** * * Read the control keywords. * ************************/ nfields = sscanf( inputline, "%s %s %s %s %s %s", keyword, param1, param2, param3, param4, param5); if( nfields < 1 ) strcpy( keyword, ""); if( strcmp( keyword, "END") == 0) break; else if( strcmp( keyword, "") == 0) x = 1.0; /* blank line */ else if( strcmp( keyword, "COMMENT=") == 0) x = 1.0; /* comment line */ else if( strcmp( keyword, "STAR1=") == 0) strcpy( mode.star1, param1); else if( strcmp(keyword, "STAR1SPOTS=") == 0) strcpy( mode.star1spots, param1); else if( strcmp( keyword, "STAR2=") == 0) strcpy( mode.star2, param1); else if( strcmp(keyword, "STAR2SPOTS=") == 0) strcpy( mode.star2spots, param1); else if( strcmp( keyword, "DISK=") == 0) strcpy( mode.disk, param1); else if( strcmp( keyword, "THIRDLIGHT=") == 0) strcpy( mode.thirdlight, param1); else if( strcmp( keyword, "REFLECTION=") == 0) strcpy( mode.reflection, param1); else if( strcmp( keyword, "WRITEGRIDS=") == 0) strcpy( mode.writegrids, param1); else if( strcmp( keyword, "RUNTYPE=") == 0) { strcpy( mode.runtype, param1); strcpy( mode.plotprogram, param2); } else if( strcmp( keyword, "NORMALIZE=") == 0) { strcpy( mode.normalize, param1); if( strcmp( mode.normalize, "SETVALUE") == 0) { if( nfields < 4 ) quit("Too few parameters for keyword NORMALIZE."); sscanf( param2, "%lf", &starsys.normphase); sscanf( param3, "%lf", &starsys.normflux); } if( strcmp( mode.normalize, "MAXVALUE") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword NORMALIZE."); sscanf( param2, "%lf", &starsys.normflux); } if( strcmp( mode.normalize, "LSQR") == 0) { if( nfields < 4 ) quit("Too few parameters for keyword NORMALIZE."); sscanf( param2, "%lf", &mode.fitphase[0] ); sscanf( param3, "%lf", &mode.fitphase[1] ); } } else if( strcmp( keyword, "PHASESOURCE=") == 0) { if( strcmp( param1, "EVENSPACE") == 0) { if( nfields < 5 ) quit("Too few parameters for keyword PHASESOURCE."); sscanf( param2, "%lf", &firstphase); sscanf( param3, "%lf", &lastphase); sscanf( param4, "%lf", &deltaphase); OrbPhases( param1, firstphase, lastphase, deltaphase); } else if( strcmp( param1, "FILE") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword PHASESOURCE."); OrbPhases( param2, firstphase, lastphase, deltaphase); } else quit("PHASESOURCE must be EVENSPACE or FILE."); } else if( strcmp( keyword, "READDATA=") == 0) { strcpy( mode.readdata, param1 ); if( strcmp( mode.readdata, "ON" ) == 0) { if( nfields < 3 ) quit("Too few parameters for keyword READDATA."); strcpy( mode.datafile, param2 ); } } /**************************** * * Read the keywords referring to the whole system. * ******************************/ else if( strcmp( keyword, "Q=") == 0) sscanf( param1, "%lf", &starsys.q); else if( strcmp( keyword, "MASS1=") == 0) sscanf( param1, "%lf", &starsys.mass1); else if( strcmp( keyword, "INCLINATION=") == 0) sscanf( param1, "%lf", &starsys.inclination); else if( strcmp( keyword, "FILTERWAVEL=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword FILTERWAVEL."); strcpy( starsys.filter, param1); sscanf( param2, "%lf", &starsys.wavelength); } else if( strcmp( keyword, "PHASEOFFSET=") == 0) sscanf( param1, "%lf", &starsys.phaseoffset); /********************* * * Read keywords referring to star 1. * **********************/ else if( strcmp( keyword, "THETAPOINTS1=") == 0) sscanf( param1, "%ld", &star1.thetapoints); else if( strcmp( keyword, "PHIPOINTS1=") == 0) sscanf( param1, "%ld", &star1.phipoints); else if( strcmp( keyword, "FRACRADIUS1=") == 0) sscanf( param1, "%lf", &star1.fracRadius); else if( strcmp( keyword, "GRAVBETA1=") == 0) sscanf( param1, "%lf", &star1.gravbeta); else if( strcmp( keyword, "POLETEMP1=") == 0) sscanf( param1, "%lf", &star1.poletemp); else if( strcmp( keyword, "FLUXSOURCE1=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword FLUXSOURCE1."); strcpy( star1.fluxsource, param1); if( strcmp( star1.fluxsource, "BB") == 0) sscanf( param2, "%lf", &star1.bblimbdark); else if( strcmp( star1.fluxsource, "ITABLE") == 0) sscanf( param2, "%lf", &star1.logg); else quit("Unrecognized flux source for star 1."); } else if( strcmp( keyword, "ALBEDO1=") == 0) sscanf( param1, "%lf", &star1.albedo); /*************************************************** * * Keywords concerning star 1 spots. * *****************************************************/ else if( strcmp(keyword, "STAR1SPOT=") == 0) { if( nfields < 5 ) quit("Too few parameters in keyword STAR1SPOT."); star1spot.nspots += 1; if( star1spot.nspots >= 20 ) quit("Too many star 1 spots."); sscanf( param1, "%lf", &star1spot.theta[star1spot.nspots]); sscanf( param2, "%lf", &star1spot.phi[star1spot.nspots]); sscanf( param3, "%lf", &star1spot.radius[star1spot.nspots]); sscanf( param4, "%lf", &star1spot.SpotToverStarT[star1spot.nspots]); star1spot.theta[star1spot.nspots] *= TWOPI / 360.0; star1spot.phi[star1spot.nspots] *= TWOPI / 360.0; star1spot.radius[star1spot.nspots] *= TWOPI / 360.0; } /******************* * * Read keywords referring to star 2. * *********************/ else if( strcmp( keyword, "THETAPOINTS2=") == 0) sscanf( param1, "%ld", &star2.thetapoints); else if( strcmp( keyword, "PHIPOINTS2=") == 0) sscanf( param1, "%ld", &star2.phipoints); else if( strcmp( keyword, "FRACRADIUS2=") == 0) sscanf( param1, "%lf", &star2.fracRadius); else if( strcmp( keyword, "GRAVBETA2=") == 0) sscanf( param1, "%lf", &star2.gravbeta); else if( strcmp( keyword, "POLETEMP2=") == 0) sscanf( param1, "%lf", &star2.poletemp); else if( strcmp( keyword, "FLUXSOURCE2=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword FLUXSOURCE2."); strcpy( star2.fluxsource, param1); if( strcmp( star2.fluxsource, "BB") == 0) sscanf( param2, "%lf", &star2.bblimbdark); else if( strcmp( star2.fluxsource, "ITABLE") == 0) sscanf( param2, "%lf", &star2.logg); else quit("Unrecognized flux source for star 2."); } else if( strcmp( keyword, "ALBEDO2=") == 0) sscanf( param1, "%lf", &star2.albedo); /*************************************************** * * Keywords concerning star 2 spots. * *****************************************************/ else if( strcmp(keyword, "STAR2SPOT=") == 0) { if( nfields < 5 ) quit("Too few parameters in keyword STAR2SPOT."); star2spot.nspots += 1; if( star2spot.nspots >= 20 ) quit("Too many star 2 spots."); sscanf( param1, "%lf", &star2spot.theta[star2spot.nspots]); sscanf( param2, "%lf", &star2spot.phi[star2spot.nspots]); sscanf( param3, "%lf", &star2spot.radius[star2spot.nspots]); sscanf( param4, "%lf", &star2spot.SpotToverStarT[star2spot.nspots]); star2spot.theta[star2spot.nspots] *= TWOPI / 360.0; star2spot.phi[star2spot.nspots] *= TWOPI / 360.0; star2spot.radius[star2spot.nspots] *= TWOPI / 360.0; } /*********************** * * Read keywords referring to the entire disk. * ************************/ else if( strcmp( keyword, "RPOINTSDISK=") == 0) { sscanf( param1, "%ld", &disktop.rpoints); sscanf( param1, "%ld", &diskrim.rpoints); } else if( strcmp( keyword, "PHIPOINTSDISK=") == 0) { sscanf( param1, "%ld", &disktop.phipoints); sscanf( param1, "%ld", &diskrim.phipoints); } else if( strcmp( keyword, "RDISKIN=") == 0) sscanf( param1, "%lf", &disktop.innerR); else if( strcmp( keyword, "RDISKOUT=") == 0) { sscanf( param1, "%lf", &disktop.outerR); sscanf( param1, "%lf", &diskrim.outerR); } else if( strcmp( keyword, "FLAREANGLE=") == 0) sscanf( param1, "%lf", &disktop.flareangle); else if( strcmp( keyword, "DISKTRANSMIT=") == 0) sscanf( param1, "%lf", &disktop.transmission); /*********************** * * Read keywords referring to the top of the disk. * ************************/ else if( strcmp( keyword, "TOPTEMP=") == 0) { strcpy( disktop.Tsource, param1); if( strcmp( disktop.Tsource, "POLYNOMIAL") == 0) { if( nfields < 6 ) quit("Too few parameters for keyword TEMPDISK."); sscanf( param2, "%lf", &disktop.Tcoef[0] ); sscanf( param3, "%lf", &disktop.Tcoef[1] ); sscanf( param4, "%lf", &disktop.Tcoef[2] ); sscanf( param5, "%lf", &disktop.Tcoef[3] ); } else if( strcmp( disktop.Tsource, "STEADYSTATE") == 0) { if( nfields < 4 ) quit("Too few parameters for keyword TEMPDISK."); sscanf( param2, "%lf", &disktop.Tcoef[0] ); sscanf( param3, "%lf", &disktop.Tcoef[1] ); } else if( strcmp( disktop.Tsource, "FILE") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword TEMPDISK."); strcpy( disktop.Tfile, param2); } else { quit("Invalid temperature source for TEMPDISK."); } } else if( strcmp( keyword, "TOPFLUXSOURCE=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword TOPFLUXSOURCE."); strcpy( disktop.fluxsource, param1); if( strcmp( disktop.fluxsource, "BB") == 0) sscanf( param2, "%lf", &disktop.bblimbdark); else if( strcmp(disktop.fluxsource, "ITABLE") == 0) sscanf( param2, "%lf", &disktop.logg); else quit("Unrecognized flux source for the disk top."); } else if( strcmp( keyword, "TOPSPOTTEMP=") == 0) { sscanf( param1, "%lf", &disktop.spotT ); } else if( strcmp( keyword, "TOPSPOTR=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword TOPSPOTR."); sscanf( param1, "%lf", &disktop.spotr[0] ); sscanf( param2, "%lf", &disktop.spotr[1] ); } else if( strcmp( keyword, "TOPSPOTPHI=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword TOPSPOTPHI."); sscanf( param1, "%lf", &disktop.spotphi[0] ); sscanf( param2, "%lf", &disktop.spotphi[1] ); } /*********************** * * Read keywords referring to the rim of the disk. * ************************/ else if( strcmp( keyword, "RIMTEMP=") == 0) { sscanf( param1, "%lf", &diskrim.T ); } else if( strcmp( keyword, "RIMFLUXSOURCE=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword RIMFLUXSOURCE."); strcpy( diskrim.fluxsource, param1); if( strcmp( diskrim.fluxsource, "BB") == 0) sscanf( param2, "%lf", &diskrim.bblimbdark); else if( strcmp(diskrim.fluxsource, "ITABLE") == 0) sscanf( param2, "%lf", &diskrim.logg); else quit("Unrecognized flux source for the disk rim."); } else if( strcmp( keyword, "RIMSPOTTEMP=") == 0) { sscanf( param1, "%lf", &diskrim.spotT ); } else if( strcmp( keyword, "RIMSPOTPHI=") == 0) { if( nfields < 3 ) quit("Too few parameters for keyword RIMSPOTPHI."); sscanf( param1, "%lf", &diskrim.spotphi[0] ); sscanf( param2, "%lf", &diskrim.spotphi[1] ); } /*********************** * * Read keywords referring to the "third light". * ************************/ else if( strcmp( keyword, "FRACTION3=") == 0) { sscanf( param1, "%lf", &thirdlight.fraction ); } else if( strcmp( keyword, "PHASE3=") == 0) { sscanf( param1, "%lf", &thirdlight.orbphase ); } /********************** * * Unrecognized keyword. * ***********************/ else { printf(" Unrecognized keyword in readParams!\n"); printf(" keyword =%20s\n", keyword); quit(""); } } fclose( in ); return; } void OrbPhases( char *phasesource, double firstphase, double lastphase, double deltaphase) /************************* * * This function makes the array data.phases[], which contains the * orbital phases at which the light curve will be calculated. * Input Data: * phasesource = either the word "EVENSPACE" or "filename" * firstphase, * lastphase, * deltaphase = obvious meanings, but only used if phasesource * is "EVENSPACE" * * NOTE: all phases from -0.5 to 1.0 are okay, but phases * between 0.5 and 1.0 will be changed to equivalent phases * between -0.5 and 0.5 in the main program. * * if phasesource = FILE * the phases will be read from the file whose name is in * mode.phasefile[], at one orbital phase per line in the file. * * NOTE: the phase is assumed to be the first number in the line, * and all other information on the line is discarded. This means * that the phases can be read from the data file if desired. * * Output Data: * calc.npoints * calc.phases[] * **************************/ { FILE *in; char inputline[81]; long int i, j; double x; if( strcmp( phasesource, "EVENSPACE") == 0) { if( deltaphase == 0.0 ) quit("deltaphase cannot be 0.0 if PHASESOURCE= EVENSPACE"); if( lastphase < firstphase ) quit("lastphase must be greater than or equal to firstphase."); for( j = 0; ; j++){ if( j >= PHASEPOINTS ) quit("Too many orbital phases in synthetic light curve."); x = firstphase + j * deltaphase; if( x > lastphase ) break; if( (x < -0.5) || (x > 1.0) ) quit("Orbital phases out of range."); calc.phases[j] = x; calc.npoints = j + 1; } } else { if( ( in = fopen( phasesource, "r")) == NULL) { printf(" Cannot open file %s which is supposed\n", phasesource); quit("to contain orbital phases."); } calc.npoints = 0; while( fgets( inputline, 80, in ) != NULL ) { sscanf( inputline, "%lf", &x); if( calc.npoints > PHASEPOINTS ) quit("Too many orbital phases in synthetic light curve."); if( (x < -0.5) || (x > 1.0) ) { printf("Orbital phase out of range in file %s\n.", phasesource); quit(""); } calc.phases[calc.npoints] = x; calc.npoints = calc.npoints + 1; } fclose( in ); } /********************************* * * Make sure the phases are in order of increasing phase. * **********************************/ if( calc.npoints <= 1 ) return; for( i = 0; i < calc.npoints; i++) { for( j = i+1; j < calc.npoints; j++) { if( calc.phases[j] < calc.phases[i] ) { x = calc.phases[i]; calc.phases[i] = calc.phases[j]; calc.phases[j] = x; } } } return; } void readData( char *filename, struct lightcurve *data ) /***************************** * * This function reads a data lightcurve from a file. * Input Data: * *filename = pointer to a char array containing the file name * The file format should have one data point per * line, each line with three floating point numbers * arranged thus: * phase flux stand-dev-of-flux * * The standard deviations must all be present or all * be missing. If some are present and some missing * the function issues a message and halts. * * *data = pointer to a lightcurve structure into which the * data light curve is to be put. * * Output data: * (*data).npoints = number of data points * (*data).phases[] = the phases of the points * (*data).fluxes[] = the fluxes * (*data).stdev[] = the standard deviations. * ******************************/ { FILE *in; char inputline[81]; long int npoints, i, j, nfields, minfields, maxfields; double phase, flux, stdev, x; if( (in = fopen( filename, "r" )) == NULL ) { printf(" Cannot open file %s which is supposed\n", filename); quit("to contain the observed light curve."); } npoints = 0; minfields = 10; maxfields = 0; while( fgets( inputline, 80, in) != NULL ) { npoints = npoints + 1; if( npoints > PHASEPOINTS ) quit("Too many phase points in observed light curve."); nfields = sscanf( inputline, "%lf %lf %lf", &phase, &flux, &stdev ); if( nfields < minfields ) minfields = nfields; if( nfields > maxfields ) maxfields = nfields; if( (phase < -0.5) || (phase > 1.0) ) { printf(" Phase of data point %3ld out of range.\n", npoints); quit(""); } (*data).phases[npoints - 1] = phase; (*data).fluxes[npoints - 1] = flux; if( nfields >= 3 ) { if( stdev <= 0.0 ) { printf(" Stand. dev. of data point %3ld out of range.\n", npoints); quit(""); } (*data).stdev[npoints - 1] = stdev; } else { (*data).stdev[npoints - 1] = 1.0; } } if( npoints <= 0 ) quit("No data points in the file containing the observed light curve."); (*data).npoints = npoints; if( maxfields != minfields ) quit("One or more standard deviations missing in data file."); /******************* * * Put the points in the observed light curve in order of * increasing orbital phase. * ********************/ if( npoints == 1 ) return; for( i = 0; i < (npoints-1); i++) { for( j = i+1; j < npoints; j++) { if( (*data).phases[j] < (*data).phases[i] ) { x = (*data).phases[i]; (*data).phases[i] = (*data).phases[j]; (*data).phases[j] = x; x = (*data).fluxes[i]; (*data).fluxes[i] = (*data).fluxes[j]; (*data).fluxes[j] = x; x = (*data).stdev[i]; (*data).stdev[i] = (*data).stdev[j]; (*data).stdev[j] = x; } } } return; } long int ReadItable( char *filter, double logg, double *T, double *Io, double *a, double *b) /************************************** * * This function reads the intensity table and extracts only those parts * with the correct filter and log g. The entries are then sorted in * order of increasing temperature. * * It returns the number of entries found. * **************************************/ { char *filename, inputline[81], inputfilter[20]; double temp, inputlogg; FILE *in; long int i, j, n; filename = "itable.dat"; if( ( in = fopen(filename, "r") ) == NULL) { printf(" Cannot open ITABLE file %s\n", filename); quit(""); } n = 0; while( fgets( inputline, 80, in ) != NULL ) { sscanf( inputline, "%s %lf %lf %lf %lf %lf", inputfilter, &inputlogg, &T[n], &Io[n], &a[n], &b[n] ); if( strcmp( inputfilter, "*") == 0 ) { n = n; /* skip over the comment line. */ } else if( strcmp( inputfilter, filter) == 0 ) { if( fabs(inputlogg - logg) <= 0.1 ) n = n + 1; } } fclose( in); if( n == 0 ) { printf("\n No intensities found for filter %s and log g = %4.2lf\n", filter, logg); quit(""); } if( n > 1 ) { for( i = 0; i < (n-1); i++ ) { for( j = i; j < n; j++ ) { if( T[j] < T[i] ) { temp = T[j]; T[j] = T[i]; T[i] = temp; temp = Io[j]; Io[j] = Io[i]; Io[i] = temp; temp = a[j]; a[j] = a[i]; a[i] = temp; temp = b[j]; b[j] = b[i]; b[i] = temp; } } } } return( n ); } void readDiskT( char *filename ) /************************* * * This function reads disk temperatures from a file. * In this version the disk temperature is assumed to be a function * only of R (not phi) and is given in a file containing a two-column * ascii table: * r1 t1 * r2 t2 * . . * . . * . . * The radii are in units of the distance from the center of star 1 * to the inner Lagrangian point, r/RL1, and are the radii in the * plane of the disk (not along the surface of the disk. * * The function linearly interpolates in the table to get temperatures * at the grid points on the disk. Grid points at radii outside * the bounds of the table are set equal to the values at the limits * of the table. * *************************/ { FILE *in; double r[LINEPOINTS], t[LINEPOINTS], slope, massfraction, L1distance, twopi, projectfactor, projectedr; long int tpoints, ir, irmax, j, iphi, iphimax; if( (in = fopen( filename, "r" )) == NULL ) quit("Cannot open file with disk T(R) data."); for( tpoints = 0; tpoints < LINEPOINTS; tpoints++) { if( fscanf(in, "%lf %lf", &r[tpoints], &t[tpoints]) == EOF) break; } fclose( in ); /************************ * * Put r[] in units of the separation of the two stars. * ************************/ massfraction = starsys.q / ( 1.0 + starsys.q ); L1distance = findL1( massfraction ); for( j = 0; j < tpoints; j++) { r[j] = r[j] * L1distance; } twopi = 2.0 * 3.14159265358979; projectfactor = cos( twopi * (disktop.flareangle / 360.0 ) ); irmax = disktop.rpoints - 4; iphimax = disktop.phipoints - 1; for( ir = 0; ir <= irmax; ir++) { projectedr = projectfactor * rdisk[ir][0]; if( projectedr < r[0] ) { localTdisk[ir][0] = t[0]; } else if( projectedr >= r[tpoints-1] ) { localTdisk[ir][0] = t[tpoints-1]; } else { for( j = 1; j < tpoints; j++) { if( (projectedr >= r[j-1]) && (projectedr < r[j]) ) { slope = ( t[j] - t[j-1] ) / ( r[j] - r[j-1] ); localTdisk[ir][0] = slope * ( projectedr - r[j-1] ) + t[j-1]; } } } for( iphi = 1; iphi <= iphimax; iphi++) { localTdisk[ir][iphi] = localTdisk[ir][0]; } } return; }