/* 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 <math.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#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);
}
syntax highlighted by Code2HTML, v. 0.9.1