/* NIGHTFALL Light Curve Synthesis Program */ /* Copyright (C) 1998 Rainer Wichmann */ /* */ /* This program is free software; you can redistribute it */ /* and/or modify */ /* it under the terms of the GNU General Public License as */ /* published by */ /* the Free Software Foundation; either version 2 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program; if not, write to the Free Software */ /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include #include #include #include "Light.h" /****************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Define parameters of binary geometry (non-overcontact) @param (int) n The stellar component @return (int) The exit status @heading Star Setup *******************************************************************/ int DefineParam(int n) { double RadPower; /* for volume computation */ double RadPol; /* polar radius */ double NN; /* (Binary[n].Mq + 1.0) / 2.0 */ int test = 0; /* exit status of RootFind */ BinaryComponent *BinPtr; /* pointer to binary */ BinPtr = &Binary[n]; Mq = BinPtr->Mq; /* set mass ratio to appropriate value */ F = 1.0; /* set nonsync rot to 1 for Lagrange Points */ /* ----------- Lagrange one (on x-axis) ----------------------- */ BinPtr->RLag1 = RootFind(LagrangeOne, (FLT_EPSILON), (1.0 - FLT_EPSILON), 1.0e-7, "DefineParam", &test); if (test == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) { fprintf(stderr, "\n"); fprintf(stderr, _(" Star: %d \n"), n+1); fprintf(stderr, _("LagrangeOne is at %8.6f\n"), BinPtr->RLag1); } if (BinPtr->Mq <= 1.0) { /* -------------------- Lagrange Two ---------------------- */ BinPtr->RLag2 = RootFind(LagrangeTwo, (1.0 + FLT_EPSILON), 2.0, 1.0e-7, "DefineParam", &test); if (test == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeTwo is at %8.6f\n"), BinPtr->RLag2); /* ----------------- Lagrange Three -------------------------- */ BinPtr->RLag3 = RootFind(LagrangeThree, (0.0 - FLT_EPSILON), -2.0, 1.0e-7, "DefineParam", &test); if (test == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeThree is at %8.6f\n"), BinPtr->RLag3); } else { /* ------------------ Lagrange Two -------------------------- */ BinPtr->RLag2 = RootFind(LagrangeThree, (0.0 - FLT_EPSILON), -2.0, 1.0e-7, "DefineParam", &test); if (test == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeTwo is at %8.6f\n"), BinPtr->RLag2); /* -------------------- Lagrange Three ----------------------- */ BinPtr->RLag3 = RootFind(LagrangeTwo, (1.0 + FLT_EPSILON), 2.0, 1.0e-7, "DefineParam", &test); if (test == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeThree is at %8.6f\n"), BinPtr->RLag3); } /* -------------- Potential at Lagrange Points -------------- */ /* F is global variable */ if (n == Primary && Flags.asynchron1 == ON) F = BinPtr->Fratio; /* set nonsync rot to appropriate value */ else if (n == Secondary && Flags.asynchron2 == ON) F = BinPtr->Fratio; else F = 1.0; BinPtr->RCLag1 = 1.0 / BinPtr->RLag1 + BinPtr->Mq*(1.0/(1.0-BinPtr->RLag1)-BinPtr->RLag1) + F*F*(BinPtr->RLag1*BinPtr->RLag1)*(BinPtr->Mq+1.0)/2.0; if (BinPtr->Mq <= 1.0) { BinPtr->RCLag2 = 1.0 / BinPtr->RLag2 + BinPtr->Mq * (1.0/(BinPtr->RLag2-1.0)-BinPtr->RLag2) + F*F*(BinPtr->RLag2*BinPtr->RLag2)*(BinPtr->Mq+1.0)/2.0; BinPtr->RCLag3 = 1.0 / fabs(BinPtr->RLag3) + BinPtr->Mq * (1.0/(1.0+fabs(BinPtr->RLag3)) - BinPtr->RLag3) + F*F*(BinPtr->RLag3*BinPtr->RLag3)*(BinPtr->Mq+1.0)/2.0; } else { BinPtr->RCLag2 = 1.0 / fabs(BinPtr->RLag2) + BinPtr->Mq * (1.0/(1.0+fabs(BinPtr->RLag2)) - BinPtr->RLag2) + F*F*(BinPtr->RLag2*BinPtr->RLag2)*(BinPtr->Mq+1.0)/2.0; BinPtr->RCLag3 = 1.0 / BinPtr->RLag3 + BinPtr->Mq * (1.0/(BinPtr->RLag3-1.0)-BinPtr->RLag3) + F*F*(BinPtr->RLag3*BinPtr->RLag3)*(BinPtr->Mq+1.0)/2.0; } /* Critical Radius on X axis, can be different from */ /* R(LagrangeOne) in case of nonsynchroneous rotation */ if ( (Flags.asynchron1 == ON && n == Primary) || (Flags.asynchron2 == ON && n == Secondary) ){ /* not used if overcontact system */ F = BinPtr->Fratio; BinPtr->RXCrit = RootFind(LagrangeOne, 0.00001, 0.99999, 1.0e-7, "DefineParam", &test); if (test == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Critical Radius on X-Axis is %8.6f\n"), BinPtr->RXCrit); } else { BinPtr->RXCrit = BinPtr->RLag1; } /* --- Roche Potential for max lobe, ie. on x-axis towards L1 */ BinPtr->RochePotCrit = 1.0 / BinPtr->RXCrit + BinPtr->Mq*(1.0/(1.0-BinPtr->RXCrit)-BinPtr->RXCrit) + F*F*(BinPtr->RXCrit*BinPtr->RXCrit) *(BinPtr->Mq+1.0)/2.0; if ((Flags.debug[GEO] == ON)) fprintf(stderr, _("Critical Potential is %8.6f\n"), BinPtr->RochePotCrit); /* --- RCrit (on z axis) , ie. critical Polar Radius ------- */ RochePot = BinPtr->RochePotCrit; lambda = 0.0; nu = 1.0; BinPtr->RCrit = RootFind(RocheSurface, 0.00001, 0.99999, 1.0e-7, "DefineParam", &test); if (test == 1) return(8); /* ------------ DEFINE ACTUAL POLAR RADIUS ------------------- */ /* AND CORRESPONDING POTENTIAL */ BinPtr->Radius = BinPtr->RocheFill * BinPtr->RCrit; if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Critical Polar Radius: %8.6f, Polar Radius: %16.12f\n"), BinPtr->RCrit, BinPtr->Radius); BinPtr->RochePot = 1.0/BinPtr->Radius + BinPtr->Mq/sqrt(1 + (BinPtr->Radius*BinPtr->Radius)); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Roche Potential: %8.6f\n"), BinPtr->RochePot); /* mean radius as defined by Kopal */ BinPtr->RadMean = 1.0 / (BinPtr->RochePot - BinPtr->Mq); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Mean Radius: %8.6f\n"), BinPtr->RadMean); /* Volume analytic as by Kopal exact to O(r**11) */ NN = (BinPtr->Mq + 1.0) / 2.0; RadPower = BinPtr->RadMean * SQR(BinPtr->RadMean); /* r^3 */ BinPtr->Volume = 1.0 + 2.0 * NN * RadPower; RadPower = SQR(RadPower); /* r^6 */ BinPtr->Volume = BinPtr->Volume + (12.0/5.0)*SQR(BinPtr->Mq)*RadPower + (32.0/5.0)*(NN*NN)*RadPower + (8.0/5.0)*NN*BinPtr->Mq*RadPower; RadPower = RadPower * SQR(BinPtr->RadMean); /* r^8 */ BinPtr->Volume = BinPtr->Volume + (15.0/7.0)*SQR(BinPtr->Mq)*RadPower; RadPower = RadPower * BinPtr->RadMean; /* r^9 */ BinPtr->Volume = BinPtr->Volume + (22.0/7.0)*SQR(BinPtr->Mq)*BinPtr->Mq*RadPower + (176.0/7.0)*(NN*NN)*NN*RadPower + (296.0/35.0)*NN*BinPtr->Mq *(2.0*BinPtr->Mq + NN)*RadPower; RadPower = RadPower * BinPtr->RadMean; /* r^10 */ BinPtr->Volume = BinPtr->Volume + (18.0/9.0)*SQR(BinPtr->Mq)*RadPower; RadPower = RadPower * BinPtr->RadMean; /* r^11 */ BinPtr->Volume = BinPtr->Volume + (157.0/7.0)*SQR(BinPtr->Mq)*BinPtr->Mq*RadPower + (26.0/35.0)*NN*BinPtr->Mq *(BinPtr->Mq + 3.0*NN)*RadPower; BinPtr->Volume = BinPtr->Volume * (4.0*PI/3.0)*SQR(BinPtr->RadMean)*BinPtr->RadMean; if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Analytic Volume: %8.6f\n"), BinPtr->Volume); /* polar radius analytic as by Kopal */ /* correct up to r**11 */ RadPol = PolRad(BinPtr->Mq, BinPtr->RadMean); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Analytic Polar Radius: %16.12f\n"), RadPol); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Analytic approximations only used in case of eccentric orbit\n")); return(0); } #ifdef HAVE_DISK /****************************************************************** @package nightfall @author @version 1.0 @short Define parameters of disk geometry (non-overcontact) @param (int) n The stellar component @return (int) The exit status @heading Disk Setup *******************************************************************/ int DefineParamDisk() { /* double RadPower;*/ /* for volume computation */ /* double RadPol; */ /* polar radius */ /* double NN; */ /* (Binary[n].Mq + 1.0) / 2.0 */ /* int test = 0; */ /* exit status of RootFind */ BinaryComponent *BinPtr; /* pointer to binary */ BinPtr = &Binary[Disk]; Mq = BinPtr->Mq; /* set mass ratio to appropriate value */ F = 1.0; /* set nonsync rot to 1 for Lagrange Points */ /* ----------- Lagrange one (on x-axis) ----------------------- */ /* if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Analytic approximations only used in case of eccentric orbit\n")); */ return(0); } #endif /****************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Define parameters of binary geometry (overcontact) @param (void) @return (int) The exit status @heading Star Setup *******************************************************************/ int DefineParamOver() { char message[128]; /* error message */ int Test = 0; /* fill factor ok ? */ int ttest = 0; /* RootFind exit status */ BinaryComponent *BinPtrP; /* pointer to Primary */ BinaryComponent *BinPtrS; /* pointer to Secondary */ double asinArg; /* argument of asin */ BinPtrP = &Binary[Primary]; BinPtrS = &Binary[Secondary]; /* --------------- Lagrange one (on x-axis) ------------------ */ Mq = BinPtrP->Mq; F = 1.0; /* set nonsync rot to 1 for Lagrange 1 */ BinPtrP->RLag1 = RootFind(LagrangeOne, FLT_EPSILON, (1.0-FLT_EPSILON), 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeOne is at %8.6f\n"), BinPtrP->RLag1); BinPtrS->RLag1 = 1.0 - BinPtrP->RLag1; BinPtrP->RXCrit = BinPtrP->RLag1; BinPtrS->RXCrit = BinPtrS->RLag1; /* --------------------- Lagrange Two --------------------- */ if (BinPtrP->Mq <= 1.0) { BinPtrP->RLag2 = RootFind(LagrangeTwo, 1.0 + FLT_EPSILON, 2.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); } else { BinPtrP->RLag2 = RootFind(LagrangeThree, 0.0 - FLT_EPSILON, -2.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); } if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeTwo is at %8.6f\n"), BinPtrP->RLag2); Mq = BinPtrS->Mq; if (BinPtrS->Mq <= 1.0) { BinPtrS->RLag2 = RootFind(LagrangeTwo, 1.0 + FLT_EPSILON, 2.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); } else { BinPtrS->RLag2 = RootFind(LagrangeThree, 0.0 - FLT_EPSILON, -2.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); } if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("LagrangeTwo (Secondary Frame) is at %8.6f\n"), BinPtrS->RLag2); Mq = BinPtrP->Mq; /* -------------------- Lagrange Three --------------------- */ if (BinPtrP->Mq <= 1.0) { BinPtrP->RLag3 = RootFind(LagrangeThree, 0.0 - FLT_EPSILON, -2.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); } else { BinPtrP->RLag3 = RootFind(LagrangeTwo, 1.0 + FLT_EPSILON, 2.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); } if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) { fprintf(stderr, _("LagrangeThree is at %8.6f\n"), BinPtrP->RLag3); fprintf(stderr, "\n"); } /* -------------- Potential at Lagrange Points -------------- */ BinPtrP->RCLag1 = 1.0 / BinPtrP->RLag1 + BinPtrP->Mq*(1.0/(1.0-BinPtrP->RLag1) -BinPtrP->RLag1) + SQR(BinPtrP->RLag1)*(BinPtrP->Mq+1.0)/2.0; if (BinPtrP->Mq <= 1.0) { BinPtrP->RCLag2 = 1.0 / BinPtrP->RLag2 + BinPtrP->Mq * (1.0/(BinPtrP->RLag2-1.0)-BinPtrP->RLag2) + SQR(BinPtrP->RLag2)*(BinPtrP->Mq+1.0)/2.0; BinPtrS->RCLag2 = 1.0 / fabs(BinPtrS->RLag2) + BinPtrS->Mq * (1.0/(1.0+fabs(BinPtrS->RLag2)) - BinPtrS->RLag2) + SQR(BinPtrS->RLag2) *(BinPtrS->Mq+1.0)/2.0; BinPtrP->RCLag3 = 1.0 / fabs(BinPtrP->RLag3) + BinPtrP->Mq * (1.0/(1.0+fabs(BinPtrP->RLag3)) - BinPtrP->RLag3) + SQR(BinPtrP->RLag3)*(BinPtrP->Mq+1.0)/2.0; } else { BinPtrP->RCLag3 = 1.0 / BinPtrP->RLag3 + BinPtrP->Mq * (1.0/(BinPtrP->RLag3-1.0)-BinPtrP->RLag3) + SQR(BinPtrP->RLag3)*(BinPtrP->Mq+1.0)/2.0; BinPtrS->RCLag2 = 1.0 / BinPtrS->RLag2 + BinPtrS->Mq * (1.0/(BinPtrS->RLag2-1.0)-BinPtrS->RLag2) + SQR(BinPtrS->RLag2) *(BinPtrS->Mq+1.0)/2.0; BinPtrP->RCLag2 = 1.0 / fabs(BinPtrP->RLag2) + BinPtrP->Mq * (1.0/(1.0+fabs(BinPtrP->RLag2)) - BinPtrP->RLag2) + SQR(BinPtrP->RLag2)*(BinPtrP->Mq+1.0)/2.0; } /* -- Roche Potential for max lobe, ie. on X-axis at L1 ------ */ BinPtrP->RochePotCrit = 1.0 / BinPtrP->RXCrit + BinPtrP->Mq * (1.0/(1.0-BinPtrP->RXCrit)-BinPtrP->RXCrit) + SQR(BinPtrP->RXCrit)*(BinPtrP->Mq+1.0)/2.0; if ((Flags.debug[GEO] == ON)) fprintf(stderr, _("Primary RochePotCrit is %8.6f\n"), BinPtrP->RochePotCrit); BinPtrS->RochePotCrit = 1.0 / BinPtrS->RXCrit + BinPtrS->Mq * (1.0/(1.0-BinPtrS->RXCrit)-BinPtrS->RXCrit) + SQR(BinPtrS->RXCrit)*(BinPtrS->Mq+1.0)/2.0; if ((Flags.debug[GEO] == ON)) fprintf(stderr, _("Secondary RochePotCrit is %8.6f\n"), BinPtrS->RochePotCrit); /* ----- RCrit (on z axis) , ie. critical Polar Radius -------- */ RochePot = BinPtrP->RochePotCrit; lambda = 0.0; nu = 1.0; Mq = BinPtrP->Mq; BinPtrP->RCrit = RootFind(RocheSurface, 0.00001, 0.999, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); /* ------------ DEFINE ACTUAL POLAR RADIUS -------------------- */ BinPtrP->Radius = BinPtrP->RocheFill*BinPtrP->RCrit; if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("RCrit: %8.6f, Polar Radius: %8.6f\n"), BinPtrP->RCrit, BinPtrP->Radius); BinPtrP->RochePot = 1.0/BinPtrP->Radius + BinPtrP->Mq/sqrt(1 + SQR(BinPtrP->Radius)); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) { fprintf(stderr, _("Primary: \n")); fprintf(stderr, _("Roche Potential At Surface: %8.6f\n"), BinPtrP->RochePot); fprintf(stderr, _("Roche Potential At L2: %8.6f\n"), BinPtrP->RCLag2); } /* if Roche Potential at surface too small, adjust fill factor */ if (BinPtrP->Mq <= 1.0) { if ((BinPtrP->RochePot - BinPtrP->RCLag2) <= FLT_EPSILON && BinPtrP->RocheFill >= 1.0002) WARNING(_("Fill factor too large, adjusting it")); while ((BinPtrP->RochePot - BinPtrP->RCLag2) <= FLT_EPSILON && BinPtrP->RocheFill >= 1.0002) { BinPtrP->RocheFill = BinPtrP->RocheFill - 0.0001; BinPtrP->Radius = BinPtrP->RocheFill*BinPtrP->RCrit; BinPtrP->RochePot = 1.0/BinPtrP->Radius + BinPtrP->Mq/sqrt(1 + SQR(BinPtrP->Radius)); if((BinPtrP->RochePot - BinPtrP->RCLag2) >= FLT_EPSILON || BinPtrP->RocheFill <= 1.0002) { sprintf(message, _("--> New Fill Factor is %7.4f"), BinPtrP->RocheFill); WARNING(message); Test = 1; } } } /* ------------ RADIUS OF SECONDARY -------------------- */ RochePot = BinPtrP->RochePot; Mq = BinPtrP->Mq; BinPtrS->Radius = RootFind(RochePerpendSecond, 0.00001, 1.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); BinPtrS->RochePot = 1.0/BinPtrS->Radius + BinPtrS->Mq/sqrt(1 + SQR(BinPtrS->Radius)); if (BinPtrP->Mq >= (1.0+DBL_EPSILON)) { if ((BinPtrS->RochePot - BinPtrS->RCLag2) <= 0.02 && BinPtrP->RocheFill >= 1.0002) WARNING(_("Fill factor too large, adjusting it")); while ((BinPtrS->RochePot - BinPtrS->RCLag2) <= 0.02 && BinPtrP->RocheFill >= 1.0002) { BinPtrP->RocheFill = BinPtrP->RocheFill - 0.0001; BinPtrP->Radius = BinPtrP->RocheFill*BinPtrP->RCrit; BinPtrP->RochePot = 1.0/BinPtrP->Radius + BinPtrP->Mq/sqrt(1 + SQR(BinPtrP->Radius)); RochePot = BinPtrP->RochePot; BinPtrS->Radius = RootFind(RochePerpendSecond, 0.00001, 1.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); BinPtrS->RochePot = 1.0/BinPtrS->Radius + BinPtrS->Mq/sqrt(1 + SQR(BinPtrS->Radius)); if((BinPtrS->RochePot - BinPtrS->RCLag2) >= 0.02 || BinPtrP->RocheFill <= 1.0002) { sprintf(message, _("--> New Fill Factor is %7.4f"), BinPtrP->RocheFill); WARNING(message); Test = 1; } } } if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Polar Radius Secondary: %8.6f\n"), BinPtrS->Radius); /* ------------ RADIUS OF THROAT -------------------- */ RochePot = BinPtrP->RochePot; Mq = BinPtrP->Mq; BinPtrP->RZatL1 = RootFind(RocheZPerpendL1, 0.00001, 1.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Radius Throat Z: %8.6f\n"), BinPtrP->RZatL1); BinPtrP->RYatL1 = RootFind(RocheYPerpendL1, 0.00001, 1.0, 1.0e-7, "DefineParamOver", &ttest); if (ttest == 1) return(8); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) fprintf(stderr, _("Radius Throat Y: %8.6f\n"), BinPtrP->RYatL1); BinPtrS->RYatL1 = BinPtrP->RYatL1; BinPtrS->RZatL1 = BinPtrP->RZatL1; /* ------------ LIMITING ANGLES AND RADII ------------------- */ BinPtrP->LimRad = sqrt( (BinPtrP->RYatL1 * BinPtrP->RYatL1) + (BinPtrP->RLag1 * BinPtrP->RLag1)); BinPtrS->LimRad = sqrt( (BinPtrS->RYatL1 * BinPtrS->RYatL1) + (BinPtrS->RLag1 * BinPtrS->RLag1)); asinArg = BinPtrP->RYatL1/BinPtrP->LimRad; BinPtrP->LimEta = asin( CLAMP(asinArg, -1.0, 1.0) ); asinArg = BinPtrS->RYatL1/BinPtrS->LimRad; BinPtrS->LimEta = asin( CLAMP(asinArg, -1.0, 1.0) ); if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) { fprintf(stderr, _("LimEta(1,2): %8.6f %8.6f, LimRad(1,2) %8.6f %8.6f,\n"), BinPtrP->LimEta,BinPtrS->LimEta, BinPtrP->LimRad,BinPtrS->LimRad); fprintf(stderr, _("Polar Radius Primary : %8.6f\n"), BinPtrP->Radius); fprintf(stderr, _("Roche Potential Primary : %8.6f\n"), BinPtrP->RochePot); fprintf(stderr, _("L2 Potential Primary : %8.6f\n"), BinPtrP->RCLag2); fprintf(stderr, _("Polar Radius Secondary: %8.6f\n"), BinPtrS->Radius); fprintf(stderr, _("Roche Potential Secondary: %8.6f\n"), BinPtrS->RochePot); fprintf(stderr, _("L2 Potential Secondary: %8.6f\n"), BinPtrS->RCLag2); } return(Test); }