1 /* NIGHTFALL Light Curve Synthesis Program */
2 /* Copyright (C) 1998 Rainer Wichmann */
3 /* */
4 /* This program is free software; you can redistribute it */
5 /* and/or modify */
6 /* it under the terms of the GNU General Public License as */
7 /* published by */
8 /* the Free Software Foundation; either version 2 of the License, or */
9 /* (at your option) any later version. */
10 /* */
11 /* This program is distributed in the hope that it will be useful, */
12 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
13 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
14 /* GNU General Public License for more details. */
15 /* */
16 /* You should have received a copy of the GNU General Public License */
17 /* along with this program; if not, write to the Free Software */
18 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
19
20
21 #include <math.h>
22 #include <stdio.h>
23 #include <string.h>
24 #include <float.h>
25 #include "Light.h"
26
27 /******************************************************************
28 @package nightfall
29 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
30 @version 1.0
31 @short Define parameters of binary geometry (non-overcontact)
32 @param (int) n The stellar component
33 @return (int) The exit status
34 @heading Star Setup
35 *******************************************************************/
DefineParam(int n)36 int DefineParam(int n)
37 {
38
39 double RadPower; /* for volume computation */
40 double RadPol; /* polar radius */
41 double NN; /* (Binary[n].Mq + 1.0) / 2.0 */
42 int test = 0; /* exit status of RootFind */
43 BinaryComponent *BinPtr; /* pointer to binary */
44
45 BinPtr = &Binary[n];
46
47 Mq = BinPtr->Mq; /* set mass ratio to appropriate value */
48 F = 1.0; /* set nonsync rot to 1 for Lagrange Points */
49
50 /* ----------- Lagrange one (on x-axis) ----------------------- */
51
52 BinPtr->RLag1 = RootFind(LagrangeOne, (FLT_EPSILON), (1.0 - FLT_EPSILON),
53 1.0e-7,
54 "DefineParam", &test);
55 if (test == 1) return(8);
56
57 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) {
58 fprintf(stderr, "\n");
59 fprintf(stderr, _(" Star: %d \n"), n+1);
60 fprintf(stderr, _("LagrangeOne is at %8.6f\n"), BinPtr->RLag1);
61 }
62
63 if (BinPtr->Mq <= 1.0) {
64
65 /* -------------------- Lagrange Two ---------------------- */
66 BinPtr->RLag2 = RootFind(LagrangeTwo, (1.0 + FLT_EPSILON), 2.0, 1.0e-7,
67 "DefineParam",
68 &test);
69 if (test == 1) return(8);
70
71 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
72 fprintf(stderr, _("LagrangeTwo is at %8.6f\n"), BinPtr->RLag2);
73
74 /* ----------------- Lagrange Three -------------------------- */
75 BinPtr->RLag3 = RootFind(LagrangeThree, (0.0 - FLT_EPSILON), -2.0,
76 1.0e-7,
77 "DefineParam", &test);
78 if (test == 1) return(8);
79
80 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
81 fprintf(stderr, _("LagrangeThree is at %8.6f\n"), BinPtr->RLag3);
82
83 } else {
84
85 /* ------------------ Lagrange Two -------------------------- */
86 BinPtr->RLag2 = RootFind(LagrangeThree, (0.0 - FLT_EPSILON), -2.0,
87 1.0e-7,
88 "DefineParam", &test);
89 if (test == 1) return(8);
90
91 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
92 fprintf(stderr, _("LagrangeTwo is at %8.6f\n"), BinPtr->RLag2);
93
94 /* -------------------- Lagrange Three ----------------------- */
95 BinPtr->RLag3 = RootFind(LagrangeTwo, (1.0 + FLT_EPSILON), 2.0,
96 1.0e-7,
97 "DefineParam", &test);
98 if (test == 1) return(8);
99
100
101 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
102 fprintf(stderr, _("LagrangeThree is at %8.6f\n"), BinPtr->RLag3);
103
104 }
105
106 /* -------------- Potential at Lagrange Points -------------- */
107
108 /* F is global variable */
109
110 if (n == Primary && Flags.asynchron1 == ON)
111 F = BinPtr->Fratio; /* set nonsync rot to appropriate value */
112 else if (n == Secondary && Flags.asynchron2 == ON)
113 F = BinPtr->Fratio;
114 else
115 F = 1.0;
116
117 BinPtr->RCLag1 = 1.0 / BinPtr->RLag1
118 + BinPtr->Mq*(1.0/(1.0-BinPtr->RLag1)-BinPtr->RLag1)
119 + F*F*(BinPtr->RLag1*BinPtr->RLag1)*(BinPtr->Mq+1.0)/2.0;
120
121 if (BinPtr->Mq <= 1.0) {
122 BinPtr->RCLag2 = 1.0 / BinPtr->RLag2
123 + BinPtr->Mq *
124 (1.0/(BinPtr->RLag2-1.0)-BinPtr->RLag2)
125 + F*F*(BinPtr->RLag2*BinPtr->RLag2)*(BinPtr->Mq+1.0)/2.0;
126 BinPtr->RCLag3 = 1.0 / fabs(BinPtr->RLag3)
127 + BinPtr->Mq *
128 (1.0/(1.0+fabs(BinPtr->RLag3)) -
129 BinPtr->RLag3)
130 + F*F*(BinPtr->RLag3*BinPtr->RLag3)*(BinPtr->Mq+1.0)/2.0;
131 } else {
132 BinPtr->RCLag2 = 1.0 / fabs(BinPtr->RLag2)
133 + BinPtr->Mq *
134 (1.0/(1.0+fabs(BinPtr->RLag2)) -
135 BinPtr->RLag2)
136 + F*F*(BinPtr->RLag2*BinPtr->RLag2)*(BinPtr->Mq+1.0)/2.0;
137 BinPtr->RCLag3 = 1.0 / BinPtr->RLag3
138 + BinPtr->Mq *
139 (1.0/(BinPtr->RLag3-1.0)-BinPtr->RLag3)
140 + F*F*(BinPtr->RLag3*BinPtr->RLag3)*(BinPtr->Mq+1.0)/2.0;
141 }
142
143
144
145 /* Critical Radius on X axis, can be different from */
146 /* R(LagrangeOne) in case of nonsynchroneous rotation */
147
148 if ( (Flags.asynchron1 == ON && n == Primary) ||
149 (Flags.asynchron2 == ON && n == Secondary) ){
150 /* not used if overcontact system */
151 F = BinPtr->Fratio;
152 BinPtr->RXCrit = RootFind(LagrangeOne, 0.00001, 0.99999, 1.0e-7,
153 "DefineParam", &test);
154 if (test == 1) return(8);
155
156 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
157 fprintf(stderr, _("Critical Radius on X-Axis is %8.6f\n"),
158 BinPtr->RXCrit);
159
160 } else {
161 BinPtr->RXCrit = BinPtr->RLag1;
162 }
163
164 /* --- Roche Potential for max lobe, ie. on x-axis towards L1 */
165
166 BinPtr->RochePotCrit = 1.0 / BinPtr->RXCrit
167 + BinPtr->Mq*(1.0/(1.0-BinPtr->RXCrit)-BinPtr->RXCrit)
168 + F*F*(BinPtr->RXCrit*BinPtr->RXCrit)
169 *(BinPtr->Mq+1.0)/2.0;
170
171 if ((Flags.debug[GEO] == ON))
172 fprintf(stderr, _("Critical Potential is %8.6f\n"), BinPtr->RochePotCrit);
173
174 /* --- RCrit (on z axis) , ie. critical Polar Radius ------- */
175
176 RochePot = BinPtr->RochePotCrit; lambda = 0.0; nu = 1.0;
177 BinPtr->RCrit = RootFind(RocheSurface, 0.00001, 0.99999, 1.0e-7,
178 "DefineParam", &test);
179 if (test == 1) return(8);
180
181 /* ------------ DEFINE ACTUAL POLAR RADIUS ------------------- */
182 /* AND CORRESPONDING POTENTIAL */
183
184 BinPtr->Radius = BinPtr->RocheFill * BinPtr->RCrit;
185 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
186 fprintf(stderr, _("Critical Polar Radius: %8.6f, Polar Radius: %16.12f\n"),
187 BinPtr->RCrit, BinPtr->Radius);
188
189 BinPtr->RochePot = 1.0/BinPtr->Radius
190 + BinPtr->Mq/sqrt(1 + (BinPtr->Radius*BinPtr->Radius));
191
192 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
193 fprintf(stderr, _("Roche Potential: %8.6f\n"), BinPtr->RochePot);
194
195 /* mean radius as defined by Kopal */
196 BinPtr->RadMean = 1.0 / (BinPtr->RochePot - BinPtr->Mq);
197 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
198 fprintf(stderr, _("Mean Radius: %8.6f\n"), BinPtr->RadMean);
199
200 /* Volume analytic as by Kopal exact to O(r**11) */
201 NN = (BinPtr->Mq + 1.0) / 2.0;
202 RadPower = BinPtr->RadMean * SQR(BinPtr->RadMean); /* r^3 */
203 BinPtr->Volume = 1.0 + 2.0 * NN * RadPower;
204 RadPower = SQR(RadPower); /* r^6 */
205 BinPtr->Volume = BinPtr->Volume
206 + (12.0/5.0)*SQR(BinPtr->Mq)*RadPower
207 + (32.0/5.0)*(NN*NN)*RadPower
208 + (8.0/5.0)*NN*BinPtr->Mq*RadPower;
209 RadPower = RadPower * SQR(BinPtr->RadMean); /* r^8 */
210 BinPtr->Volume = BinPtr->Volume
211 + (15.0/7.0)*SQR(BinPtr->Mq)*RadPower;
212 RadPower = RadPower * BinPtr->RadMean; /* r^9 */
213 BinPtr->Volume = BinPtr->Volume
214 + (22.0/7.0)*SQR(BinPtr->Mq)*BinPtr->Mq*RadPower
215 + (176.0/7.0)*(NN*NN)*NN*RadPower
216 + (296.0/35.0)*NN*BinPtr->Mq
217 *(2.0*BinPtr->Mq + NN)*RadPower;
218 RadPower = RadPower * BinPtr->RadMean; /* r^10 */
219 BinPtr->Volume = BinPtr->Volume
220 + (18.0/9.0)*SQR(BinPtr->Mq)*RadPower;
221 RadPower = RadPower * BinPtr->RadMean; /* r^11 */
222 BinPtr->Volume = BinPtr->Volume
223 + (157.0/7.0)*SQR(BinPtr->Mq)*BinPtr->Mq*RadPower
224 + (26.0/35.0)*NN*BinPtr->Mq
225 *(BinPtr->Mq + 3.0*NN)*RadPower;
226
227 BinPtr->Volume = BinPtr->Volume
228 * (4.0*PI/3.0)*SQR(BinPtr->RadMean)*BinPtr->RadMean;
229
230 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
231 fprintf(stderr, _("Analytic Volume: %8.6f\n"), BinPtr->Volume);
232
233
234 /* polar radius analytic as by Kopal */
235 /* correct up to r**11 */
236
237 RadPol = PolRad(BinPtr->Mq, BinPtr->RadMean);
238 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
239 fprintf(stderr, _("Analytic Polar Radius: %16.12f\n"), RadPol);
240
241 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
242 fprintf(stderr,
243 _("Analytic approximations only used in case of eccentric orbit\n"));
244
245 return(0);
246 }
247
248 #ifdef HAVE_DISK
249 /******************************************************************
250 @package nightfall
251 @author
252 @version 1.0
253 @short Define parameters of disk geometry (non-overcontact)
254 @param (int) n The stellar component
255 @return (int) The exit status
256 @heading Disk Setup
257 *******************************************************************/
DefineParamDisk()258 int DefineParamDisk()
259 {
260
261 /* double RadPower;*/ /* for volume computation */
262 /* double RadPol; */ /* polar radius */
263 /* double NN; */ /* (Binary[n].Mq + 1.0) / 2.0 */
264 /* int test = 0; */ /* exit status of RootFind */
265 BinaryComponent *BinPtr; /* pointer to binary */
266
267 BinPtr = &Binary[Disk];
268
269 Mq = BinPtr->Mq; /* set mass ratio to appropriate value */
270 F = 1.0; /* set nonsync rot to 1 for Lagrange Points */
271
272 /* ----------- Lagrange one (on x-axis) ----------------------- */
273 /* if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
274 fprintf(stderr,
275 _("Analytic approximations only used in case of eccentric orbit\n"));
276 */
277 return(0);
278 }
279 #endif
280
281 /******************************************************************
282 @package nightfall
283 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
284 @version 1.0
285 @short Define parameters of binary geometry (overcontact)
286 @param (void)
287 @return (int) The exit status
288 @heading Star Setup
289 *******************************************************************/
DefineParamOver()290 int DefineParamOver()
291 {
292 char message[128]; /* error message */
293 int Test = 0; /* fill factor ok ? */
294 int ttest = 0; /* RootFind exit status */
295 BinaryComponent *BinPtrP; /* pointer to Primary */
296 BinaryComponent *BinPtrS; /* pointer to Secondary */
297 double asinArg; /* argument of asin */
298
299 BinPtrP = &Binary[Primary];
300 BinPtrS = &Binary[Secondary];
301
302
303 /* --------------- Lagrange one (on x-axis) ------------------ */
304
305 Mq = BinPtrP->Mq;
306 F = 1.0; /* set nonsync rot to 1 for Lagrange 1 */
307 BinPtrP->RLag1 = RootFind(LagrangeOne, FLT_EPSILON, (1.0-FLT_EPSILON),
308 1.0e-7,
309 "DefineParamOver", &ttest);
310 if (ttest == 1) return(8);
311
312
313 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
314 fprintf(stderr, _("LagrangeOne is at %8.6f\n"), BinPtrP->RLag1);
315 BinPtrS->RLag1 = 1.0 - BinPtrP->RLag1;
316
317 BinPtrP->RXCrit = BinPtrP->RLag1;
318 BinPtrS->RXCrit = BinPtrS->RLag1;
319
320 /* --------------------- Lagrange Two --------------------- */
321
322 if (BinPtrP->Mq <= 1.0) {
323 BinPtrP->RLag2 = RootFind(LagrangeTwo, 1.0 + FLT_EPSILON, 2.0, 1.0e-7,
324 "DefineParamOver", &ttest);
325 if (ttest == 1) return(8);
326 } else {
327 BinPtrP->RLag2 = RootFind(LagrangeThree, 0.0 - FLT_EPSILON, -2.0, 1.0e-7,
328 "DefineParamOver", &ttest);
329 if (ttest == 1) return(8);
330 }
331
332 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
333 fprintf(stderr, _("LagrangeTwo is at %8.6f\n"), BinPtrP->RLag2);
334
335 Mq = BinPtrS->Mq;
336
337 if (BinPtrS->Mq <= 1.0) {
338 BinPtrS->RLag2 = RootFind(LagrangeTwo, 1.0 + FLT_EPSILON, 2.0, 1.0e-7,
339 "DefineParamOver", &ttest);
340 if (ttest == 1) return(8);
341 } else {
342 BinPtrS->RLag2 = RootFind(LagrangeThree, 0.0 - FLT_EPSILON, -2.0, 1.0e-7,
343 "DefineParamOver", &ttest);
344 if (ttest == 1) return(8);
345 }
346
347 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
348 fprintf(stderr, _("LagrangeTwo (Secondary Frame) is at %8.6f\n"),
349 BinPtrS->RLag2);
350
351 Mq = BinPtrP->Mq;
352
353 /* -------------------- Lagrange Three --------------------- */
354
355 if (BinPtrP->Mq <= 1.0) {
356 BinPtrP->RLag3 = RootFind(LagrangeThree, 0.0 - FLT_EPSILON, -2.0, 1.0e-7,
357 "DefineParamOver", &ttest);
358 if (ttest == 1) return(8);
359 } else {
360 BinPtrP->RLag3 = RootFind(LagrangeTwo, 1.0 + FLT_EPSILON, 2.0, 1.0e-7,
361 "DefineParamOver", &ttest);
362 if (ttest == 1) return(8);
363 }
364
365 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) {
366 fprintf(stderr, _("LagrangeThree is at %8.6f\n"), BinPtrP->RLag3);
367 fprintf(stderr, "\n");
368 }
369
370 /* -------------- Potential at Lagrange Points -------------- */
371
372 BinPtrP->RCLag1 = 1.0 / BinPtrP->RLag1
373 + BinPtrP->Mq*(1.0/(1.0-BinPtrP->RLag1)
374 -BinPtrP->RLag1)
375 + SQR(BinPtrP->RLag1)*(BinPtrP->Mq+1.0)/2.0;
376
377 if (BinPtrP->Mq <= 1.0) {
378 BinPtrP->RCLag2 = 1.0 / BinPtrP->RLag2
379 + BinPtrP->Mq *
380 (1.0/(BinPtrP->RLag2-1.0)-BinPtrP->RLag2)
381 + SQR(BinPtrP->RLag2)*(BinPtrP->Mq+1.0)/2.0;
382
383 BinPtrS->RCLag2 = 1.0 / fabs(BinPtrS->RLag2)
384 + BinPtrS->Mq *
385 (1.0/(1.0+fabs(BinPtrS->RLag2))
386 - BinPtrS->RLag2)
387 + SQR(BinPtrS->RLag2)
388 *(BinPtrS->Mq+1.0)/2.0;
389
390 BinPtrP->RCLag3 = 1.0 / fabs(BinPtrP->RLag3)
391 + BinPtrP->Mq *
392 (1.0/(1.0+fabs(BinPtrP->RLag3)) -
393 BinPtrP->RLag3)
394 + SQR(BinPtrP->RLag3)*(BinPtrP->Mq+1.0)/2.0;
395 } else {
396 BinPtrP->RCLag3 = 1.0 / BinPtrP->RLag3
397 + BinPtrP->Mq *
398 (1.0/(BinPtrP->RLag3-1.0)-BinPtrP->RLag3)
399 + SQR(BinPtrP->RLag3)*(BinPtrP->Mq+1.0)/2.0;
400
401 BinPtrS->RCLag2 = 1.0 / BinPtrS->RLag2
402 + BinPtrS->Mq *
403 (1.0/(BinPtrS->RLag2-1.0)-BinPtrS->RLag2)
404 + SQR(BinPtrS->RLag2)
405 *(BinPtrS->Mq+1.0)/2.0;
406
407 BinPtrP->RCLag2 = 1.0 / fabs(BinPtrP->RLag2)
408 + BinPtrP->Mq *
409 (1.0/(1.0+fabs(BinPtrP->RLag2)) -
410 BinPtrP->RLag2)
411 + SQR(BinPtrP->RLag2)*(BinPtrP->Mq+1.0)/2.0;
412 }
413
414 /* -- Roche Potential for max lobe, ie. on X-axis at L1 ------ */
415
416 BinPtrP->RochePotCrit = 1.0 / BinPtrP->RXCrit
417 + BinPtrP->Mq *
418 (1.0/(1.0-BinPtrP->RXCrit)-BinPtrP->RXCrit)
419 + SQR(BinPtrP->RXCrit)*(BinPtrP->Mq+1.0)/2.0;
420 if ((Flags.debug[GEO] == ON))
421 fprintf(stderr, _("Primary RochePotCrit is %8.6f\n"),
422 BinPtrP->RochePotCrit);
423
424 BinPtrS->RochePotCrit = 1.0 / BinPtrS->RXCrit
425 + BinPtrS->Mq *
426 (1.0/(1.0-BinPtrS->RXCrit)-BinPtrS->RXCrit)
427 + SQR(BinPtrS->RXCrit)*(BinPtrS->Mq+1.0)/2.0;
428 if ((Flags.debug[GEO] == ON))
429 fprintf(stderr, _("Secondary RochePotCrit is %8.6f\n"),
430 BinPtrS->RochePotCrit);
431
432 /* ----- RCrit (on z axis) , ie. critical Polar Radius -------- */
433 RochePot = BinPtrP->RochePotCrit; lambda = 0.0; nu = 1.0;
434 Mq = BinPtrP->Mq;
435 BinPtrP->RCrit = RootFind(RocheSurface, 0.00001, 0.999, 1.0e-7,
436 "DefineParamOver", &ttest);
437 if (ttest == 1) return(8);
438
439 /* ------------ DEFINE ACTUAL POLAR RADIUS -------------------- */
440
441 BinPtrP->Radius = BinPtrP->RocheFill*BinPtrP->RCrit;
442
443 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
444 fprintf(stderr, _("RCrit: %8.6f, Polar Radius: %8.6f\n"),
445 BinPtrP->RCrit, BinPtrP->Radius);
446
447 BinPtrP->RochePot = 1.0/BinPtrP->Radius
448 + BinPtrP->Mq/sqrt(1 + SQR(BinPtrP->Radius));
449
450 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) {
451 fprintf(stderr, _("Primary: \n"));
452 fprintf(stderr, _("Roche Potential At Surface: %8.6f\n"),
453 BinPtrP->RochePot);
454 fprintf(stderr, _("Roche Potential At L2: %8.6f\n"),
455 BinPtrP->RCLag2);
456 }
457
458 /* if Roche Potential at surface too small, adjust fill factor */
459
460 if (BinPtrP->Mq <= 1.0) {
461
462 if ((BinPtrP->RochePot - BinPtrP->RCLag2) <= FLT_EPSILON
463 && BinPtrP->RocheFill >= 1.0002)
464 WARNING(_("Fill factor too large, adjusting it"));
465
466 while ((BinPtrP->RochePot - BinPtrP->RCLag2) <= FLT_EPSILON
467 && BinPtrP->RocheFill >= 1.0002) {
468 BinPtrP->RocheFill = BinPtrP->RocheFill - 0.0001;
469 BinPtrP->Radius = BinPtrP->RocheFill*BinPtrP->RCrit;
470 BinPtrP->RochePot = 1.0/BinPtrP->Radius
471 + BinPtrP->Mq/sqrt(1 + SQR(BinPtrP->Radius));
472
473 if((BinPtrP->RochePot - BinPtrP->RCLag2) >= FLT_EPSILON
474 || BinPtrP->RocheFill <= 1.0002) {
475 sprintf(message, _("--> New Fill Factor is %7.4f"),
476 BinPtrP->RocheFill);
477 WARNING(message);
478 Test = 1;
479 }
480 }
481 }
482
483
484 /* ------------ RADIUS OF SECONDARY -------------------- */
485
486 RochePot = BinPtrP->RochePot;
487 Mq = BinPtrP->Mq;
488 BinPtrS->Radius = RootFind(RochePerpendSecond, 0.00001, 1.0, 1.0e-7,
489 "DefineParamOver", &ttest);
490 if (ttest == 1) return(8);
491
492 BinPtrS->RochePot = 1.0/BinPtrS->Radius
493 + BinPtrS->Mq/sqrt(1 + SQR(BinPtrS->Radius));
494
495
496 if (BinPtrP->Mq >= (1.0+DBL_EPSILON)) {
497
498 if ((BinPtrS->RochePot - BinPtrS->RCLag2) <= 0.02
499 && BinPtrP->RocheFill >= 1.0002)
500 WARNING(_("Fill factor too large, adjusting it"));
501 while ((BinPtrS->RochePot - BinPtrS->RCLag2) <= 0.02
502 && BinPtrP->RocheFill >= 1.0002) {
503 BinPtrP->RocheFill = BinPtrP->RocheFill - 0.0001;
504 BinPtrP->Radius = BinPtrP->RocheFill*BinPtrP->RCrit;
505 BinPtrP->RochePot = 1.0/BinPtrP->Radius
506 + BinPtrP->Mq/sqrt(1 + SQR(BinPtrP->Radius));
507 RochePot = BinPtrP->RochePot;
508 BinPtrS->Radius = RootFind(RochePerpendSecond,
509 0.00001, 1.0, 1.0e-7, "DefineParamOver", &ttest);
510 if (ttest == 1) return(8);
511
512 BinPtrS->RochePot = 1.0/BinPtrS->Radius
513 + BinPtrS->Mq/sqrt(1 + SQR(BinPtrS->Radius));
514 if((BinPtrS->RochePot - BinPtrS->RCLag2) >= 0.02
515 || BinPtrP->RocheFill <= 1.0002) {
516 sprintf(message, _("--> New Fill Factor is %7.4f"),
517 BinPtrP->RocheFill);
518 WARNING(message);
519 Test = 1;
520 }
521 }
522
523 }
524
525
526 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
527 fprintf(stderr, _("Polar Radius Secondary: %8.6f\n"),
528 BinPtrS->Radius);
529
530 /* ------------ RADIUS OF THROAT -------------------- */
531
532 RochePot = BinPtrP->RochePot;
533 Mq = BinPtrP->Mq;
534 BinPtrP->RZatL1 = RootFind(RocheZPerpendL1, 0.00001, 1.0, 1.0e-7,
535 "DefineParamOver", &ttest);
536 if (ttest == 1) return(8);
537
538 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
539 fprintf(stderr, _("Radius Throat Z: %8.6f\n"), BinPtrP->RZatL1);
540
541
542 BinPtrP->RYatL1 = RootFind(RocheYPerpendL1, 0.00001, 1.0, 1.0e-7,
543 "DefineParamOver", &ttest);
544 if (ttest == 1) return(8);
545
546 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON))
547 fprintf(stderr, _("Radius Throat Y: %8.6f\n"), BinPtrP->RYatL1);
548
549 BinPtrS->RYatL1 = BinPtrP->RYatL1;
550 BinPtrS->RZatL1 = BinPtrP->RZatL1;
551
552
553 /* ------------ LIMITING ANGLES AND RADII ------------------- */
554
555 BinPtrP->LimRad = sqrt( (BinPtrP->RYatL1 * BinPtrP->RYatL1)
556 + (BinPtrP->RLag1 * BinPtrP->RLag1));
557 BinPtrS->LimRad = sqrt( (BinPtrS->RYatL1 * BinPtrS->RYatL1)
558 + (BinPtrS->RLag1 * BinPtrS->RLag1));
559
560 asinArg = BinPtrP->RYatL1/BinPtrP->LimRad;
561 BinPtrP->LimEta = asin( CLAMP(asinArg, -1.0, 1.0) );
562
563 asinArg = BinPtrS->RYatL1/BinPtrS->LimRad;
564 BinPtrS->LimEta = asin( CLAMP(asinArg, -1.0, 1.0) );
565
566 if ((Flags.debug[GEO] == ON) || (Flags.debug[MINFIND] == ON)) {
567 fprintf(stderr, _("LimEta(1,2): %8.6f %8.6f, LimRad(1,2) %8.6f %8.6f,\n"),
568 BinPtrP->LimEta,BinPtrS->LimEta,
569 BinPtrP->LimRad,BinPtrS->LimRad);
570
571 fprintf(stderr,
572 _("Polar Radius Primary : %8.6f\n"), BinPtrP->Radius);
573 fprintf(stderr,
574 _("Roche Potential Primary : %8.6f\n"), BinPtrP->RochePot);
575 fprintf(stderr,
576 _("L2 Potential Primary : %8.6f\n"), BinPtrP->RCLag2);
577
578 fprintf(stderr,
579 _("Polar Radius Secondary: %8.6f\n"), BinPtrS->Radius);
580 fprintf(stderr,
581 _("Roche Potential Secondary: %8.6f\n"), BinPtrS->RochePot);
582 fprintf(stderr,
583 _("L2 Potential Secondary: %8.6f\n"), BinPtrS->RCLag2);
584 }
585 return(Test);
586 }
587
588
589
590
591
592
593
594