/******************************** * * FILE MAIN.C * ************************************/ #include "header.h" int main (void) { char yesno[20]; double phase, x, normfactor; long int i, j; /**************************** * * First set and check a few parameters internal to the program. * *****************************/ initParams(); /************************* * * Read the "starpars" file and check that the control * parameters and system parameters are valid. * *************************/ readParams(); checkModes(); checkSysParams(); checkThirdlightParams( thirdlight ); /************************ * * Read the observed light curve. * *************************/ if( strcmp( mode.readdata, "ON" ) == 0 ) readData( mode.datafile, &data1 ); /* note the pointer */ /************************** * * Calculate the surface properties of the stars and disk: * 1. Check that the parameters for the stars and disk are valid. * 2. Calculate the geometry of the stars and disk. * 3. Calculate the local effective temperatures of the stars * and disk without including irradiation. * 4. Read the necessary parts of the intensity table from a file. * 5. Calculate the local effective temperatures of the stars * including irradiation by the other star. * ****************************/ if( strcmp( mode.star1, "ON") == 0 ) { star1.massFrac = starsys.q / ( 1.0 + starsys.q ); checkStarParams( star1 ); StarGeom( &star1, r1, l1, m1, n1, gravity1, normalx1, normaly1, normalz1, dS1 ); Tstars( "star1" ); if( strcmp( star1.fluxsource, "BB") != 0 ) { Itabn1 = ReadItable( starsys.filter, star1.logg, ItabT1, ItabIo1, ItabA1, ItabB1); } } if( strcmp( mode.star1spots, "ON") == 0 ) { checkStar1SpotParams( ); for( i = 1; i <= star1spot.nspots; i++ ) { Tstar1spot( star1spot.theta[i], star1spot.phi[i], star1spot.radius[i], star1spot.SpotToverStarT[i] ); } } if( strcmp( mode.star2, "ON") == 0 ) { star2.massFrac = 1.0 / ( 1.0 + starsys.q ); checkStarParams( star2 ); StarGeom( &star2, r2, l2, m2, n2, gravity2, normalx2, normaly2, normalz2, dS2 ); Tstars( "star2" ); if( strcmp( star2.fluxsource, "BB") != 0 ) { Itabn2 = ReadItable( starsys.filter, star2.logg, ItabT2, ItabIo2, ItabA2, ItabB2); } } if( strcmp( mode.star2spots, "ON") == 0 ) { checkStar2SpotParams( ); for( i = 1; i <= star2spot.nspots; i++ ) { Tstar2spot( star2spot.theta[i], star2spot.phi[i], star2spot.radius[i], star2spot.SpotToverStarT[i] ); } } if( strcmp( mode.disk, "ON") == 0 ) { checkDiskParams(); DiskGeom(); Tdisk(); if( strcmp( disktop.fluxsource, "BB") != 0 ) { Itabntop = ReadItable( starsys.filter, disktop.logg, ItabTtop, ItabIotop, ItabAtop, ItabBtop); } if( strcmp( diskrim.fluxsource, "BB") != 0 ) { Itabnrim = ReadItable( starsys.filter, diskrim.logg, ItabTrim, ItabIorim, ItabArim, ItabBrim); } } if( strcmp( mode.reflection, "ON") == 0) irradiate(); /********************** * * Write out the file "outpars.dat", which summarizes the model, * and then, if desired, write out the surface grids for the stars * and the disk. The grids are not normally necessary. * ***********************/ writePars(); if( strcmp( mode.writegrids, "ON") == 0 ) { if( strcmp( mode.star1, "ON") == 0 ) writeStarData( star1 ); if( strcmp( mode.star2, "ON") == 0 ) writeStarData( star2 ); if( strcmp( mode.disk, "ON") == 0 ) writeDiskData(); } /************************* * * Either calculate the lightcurve * or calculate the point image of the system * or calculate the line image of the system. * **************************/ if( strcmp( mode.runtype, "LIGHTCURVE") == 0 ) { if( strcmp( mode.thirdlight, "ON" ) == 0 ) { originalflux = totalFlux( thirdlight.orbphase ); thirdflux = originalflux * ( thirdlight.fraction / (1.0 - thirdlight.fraction) ); } for( j = 0; j < calc.npoints; j++ ) { x = calc.phases[j] - starsys.phaseoffset; if( x > 0.5 ) x = x - 1.0; if( x < -0.5) x = x + 1.0; calc.fluxes[j] = totalFlux( x ) + thirdflux; /* Uncomment this line if your paranoia gets the best of you. printf(" %7.4lf\n", calc.phases[j]); */ } if( (strcmp( mode.normalize, "SETVALUE" ) == 0 ) || (strcmp( mode.normalize, "MAXVALUE" ) == 0 ) ) normfactor = Normalize( calc.npoints, calc.phases, calc.fluxes ); if( strcmp( mode.normalize, "LSQR") == 0 ) { normfactor = lsqrNorm( data1.npoints, data1.phases, data1.fluxes, data1.stdev, calc.npoints, calc.phases, calc.fluxes); starsys.chisquared = chiSqr( data1.npoints, data1.phases, data1.fluxes, data1.stdev, calc.npoints, calc.phases, calc.fluxes); writeChiSq( starsys.chisquared ); printf(" Chi^2 = %9.4le\n", starsys.chisquared); } writeLCurve( calc.npoints, calc.phases, calc.fluxes); } else if( strcmp( mode.runtype, "LINEIMAGE") == 0 ) { printf(" Binary is about to create %ld files. Is this okay (yes or no)? ", calc.npoints); scanf("%s", yesno); if( strcmp( yesno, "yes") != 0 ) quit("Good Bye."); for( j = 0; j < calc.npoints; j++ ) { x = calc.phases[j] - starsys.phaseoffset; if( x > 0.5 ) x = x - 1.0; if( x < -0.5) x = x + 1.0; lineImage( x ); /* Uncomment this line if your paranoia gets the best of you. printf(" %7.4lf\n", calc.phases[j]); */ } } else if( strcmp( mode.runtype, "POINTIMAGE") == 0 ) { printf("\n Enter the orbital phase of the image (-0.5 to 0.5): "); scanf("%lf", &phase); if( (phase < -0.5) || (phase > 0.5) ) quit("Orbital phase out of range."); pointImage( phase ); } else quit("Unrecognized RUNTYPE."); return( 0 ); } void initParams( void ) /*************************** * * Initialize all parameters. * ********************************/ { if( LINEPOINTS < 2*(SIZEDIMEN1 + SIZEDIMEN2) ) quit("LINEPOINTS in header.h not large enough."); strcpy( mode.star1, "OFF"); strcpy( mode.star2, "OFF"); strcpy( mode.star2spots, "OFF"); strcpy( mode.disk, "OFF"); strcpy( mode.thirdlight, "OFF"); strcpy( mode.reflection, "OFF"); strcpy( mode.writegrids, "OFF"); strcpy( mode.runtype, "LIGHTCURVE"); strcpy( mode.normalize, "OFF"); strcpy( mode.readdata, "OFF"); starsys.q = 0.0; starsys.mass1 = 0.0; starsys.inclination = 0.0; strcpy( starsys.filter, ""); starsys.phaseoffset = 0.0; star1.whichstar = 1; star1.thetapoints = 0; star1.phipoints = 0; star1.fracRadius = 0.0; star1.gravbeta = -1.0; star1.poletemp = 0.0; star1.albedo = 0.0; strcpy( star1.fluxsource, ""); star1spot.nspots = 0; star2.whichstar = 2; star2.thetapoints = 0; star2.phipoints = 0; star2.fracRadius = 0.0; star2.gravbeta = 0.0; star2.poletemp = 0.0; star2.albedo = 0.0; strcpy( star2.fluxsource, ""); star2spot.nspots = 0; disktop.rpoints = 0; disktop.phipoints = 0; disktop.innerR = 0.0; disktop.outerR = 0.0; disktop.flareangle = 0.0; disktop.transmission = 0.0; strcpy( disktop.Tsource, ""); disktop.Tcoef[0] = -1.0; strcpy( disktop.fluxsource, ""); disktop.bblimbdark = -1.0; disktop.logg = -1.0; disktop.spotT = -1.0; disktop.spotr[0] = -1.0; disktop.spotphi[0] = -1.0; diskrim.rpoints = 0; diskrim.phipoints = 0; diskrim.outerR = 0.0; diskrim.T = -1.0; strcpy( diskrim.fluxsource, ""); diskrim.bblimbdark = -1.0; diskrim.logg = -1.0; diskrim.spotT = -1.0; diskrim.spotphi[0] = -1.0; thirdlight.fraction = 0.0; thirdlight.orbphase = 0.0; data1.npoints = 0; calc.npoints = 0; return; }