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 <stdlib.h>
24 #include <string.h>
25 #include <float.h>
26 #include "Light.h"
27
28 /* 2002-05-03 rwichmann HIGH_PRECISION code
29 */
30 #ifdef HIGH_PRECISION
31
32 /****************************************************************************
33 @package nightfall
34 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
35 @version 1.2
36 @short Numerically differentiate dR/deta
37 @param (int) Comp The stellar component
38 @return (double) dR/deta
39 @heading Star Setup
40 ****************************************************************************/
NumDiffR(double deltaIn,double etaIn,double sinphi,double RadiusBound,double RXCrit,int * testerr,const double eps)41 double NumDiffR (double deltaIn, double etaIn, double sinphi,
42 double RadiusBound, double RXCrit, int * testerr, const double eps)
43 {
44 double val1, val0;
45 double delta;
46 double nu_store;
47 double lambda_store;
48
49 delta = eps * deltaIn;
50 delta = (delta < eps) ? eps : delta;
51
52 nu_store = nu;
53 lambda_store = lambda;
54
55 nu = sinphi * sin(etaIn-delta);
56 lambda = cos(etaIn-delta);
57 val0 = RootFind(RocheSurface, RadiusBound, RXCrit,
58 100.0 * DBL_EPSILON,
59 "NumDiffR", testerr);
60 if (*testerr == 1)
61 return(0.0);
62
63 nu = sinphi * sin(etaIn+delta);
64 lambda = cos(etaIn+delta);
65 val1 = RootFind(RocheSurface, RadiusBound, RXCrit,
66 100.0 * DBL_EPSILON,
67 "NumDiffR", testerr);
68 if (*testerr == 1)
69 return(0.0);
70
71 nu = nu_store;
72 lambda = lambda_store;
73 return ((val1-val0)/(2.0*delta));
74 }
75
76
77 /****************************************************************************
78 @package nightfall
79 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
80 @version 1.2
81 @short Numerically differentiate dR/dphi
82 @param (int) Comp The stellar component
83 @return (double) dR/dphi
84 @heading Star Setup
85 ****************************************************************************/
NumDiffRPhi(double deltaIn,double phiIn,double sineta,double RadiusBound,double RXCrit,int * testerr,const double eps)86 double NumDiffRPhi (double deltaIn, double phiIn, double sineta,
87 double RadiusBound, double RXCrit, int * testerr, const double eps)
88 {
89 double val1, val0;
90 double delta;
91 double nu_store;
92
93 delta = eps * deltaIn;
94 delta = (delta < eps) ? eps : delta;
95
96 nu_store = nu;
97
98 nu = sineta * sin(phiIn-delta);
99 val0 = RootFind(RocheSurface, RadiusBound, RXCrit,
100 100.0 * DBL_EPSILON,
101 "NumDiffRPhi", testerr);
102 if (*testerr == 1)
103 return(0.0);
104
105 nu = sineta * sin(phiIn+delta);
106 val1 = RootFind(RocheSurface, RadiusBound, RXCrit,
107 100.0 * DBL_EPSILON,
108 "NumDiffRPhi", testerr);
109 if (*testerr == 1)
110 return(0.0);
111
112 nu = nu_store;
113 return ((val1-val0)/(2.0*delta));
114 }
115 /* #ifdef HIGH_PRECISION
116 */
117 #endif
118
119 static double NormPower[NUM_COMP] = { 0.0 };
120
ResetNormPower()121 void ResetNormPower()
122 {
123 int i;
124
125 for (i = 0; i < NUM_COMP; ++i)
126 NormPower[i] = 0.0;
127
128 return;
129 }
130
131 /****************************************************************************
132 @package nightfall
133 @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
134 @version 1.2
135 @short Divide stellar surface into grid
136 @long First, set up the grid, then compute the
137 temperature distribution on the grid, including
138 spots and reflection
139 @param (int) Comp The stellar component
140 @return (int) The exit status
141 @heading Star Setup
142 ****************************************************************************/
DivideSurface(int Comp,int flag_full)143 int DivideSurface(int Comp, int flag_full)
144 {
145 register long i,j; /* loop variables */
146 long n_phi; /* # of surface steps in phi */
147 long n_eta; /* # of surface steps in eta */
148 long NumElem; /* current surface element */
149 long NextPtr; /* next surface element */
150 int testerr = 0; /* exit status of RootFind */
151
152 double eta; /* eta of surface element */
153 double phi; /* phi of surface element */
154 double sineta, sinphi; /* sin,cos of phi */
155 double coseta, cosphi; /* sin,cos of eta */
156 double rad1; /* radius of surface element */
157 double sqr_rad1; /* square of radius */
158 double delta_eta; /* step in eta */
159 double delta_phi; /* step in phi */
160 double rho; /* rho of surface element */
161 /* 2002-05-03 rwichmann HIGH_PRECISION code
162 */
163 #ifdef HIGH_PRECISION
164 double rdiff; /* dR/deta */
165 double rdiff2; /* dR/dphi */
166 double eps = pow(DBL_EPSILON, (1.0/3.0));
167 #endif
168 double cub_rho; /* cube of rho */
169 double C_eta, C_phi, C_rad1; /* surface gravity vector */
170 /* in spherical coordinates */
171 double C_x, C_y, C_z; /* surface gravity vector */
172 /* in cartesian coordinates */
173 double grav; /* surface gravity */
174 double lvec, mvec, nvec; /* surface normal vector */
175 double RXCrit; /* Critical radius */
176 double DSurf; /* surface of element */
177 double MeanGrav = 0.0; /* mean gravity */
178 double Surf = 0.0; /* total surface */
179 double MeanDarkGrav = 0.0; /* contribution to grav. darkening */
180 double RadiusBound = 0.0; /* inner bound for RootFind */
181 double T1; /* temperature of this star */
182 double T2; /* temperature of the other star */
183 double Tr; /* temperature ratio^4 */
184 double etacorr_over, eta_corr, deltaeta_corr; /* Throat adjustment */
185 long n_next, n_last, n_comp; /* Pointer setting */
186 /* for mean temperature calculation */
187 double Area; /* area of element */
188 double Weight=0.0; /* total area */
189 double TempMean=0.0; /* mean temperature */
190
191 SurfaceElement *SurfPtr; /* pointer to surface */
192 BinaryComponent *BinPtr; /* pointer to binary */
193
194
195 /* ---------- initialize ---------------------------------------------- */
196
197 SurfPtr = Surface[Comp];
198 BinPtr = &Binary[Comp];
199
200 if (flag_full == OFF)
201 goto tcalc;
202
203 Mq = Binary[Comp].Mq;
204 if (Comp == Primary && Flags.asynchron1 == ON)
205 F = Binary[Comp].Fratio;
206 else if (Comp == Secondary && Flags.asynchron2 == ON)
207 F = Binary[Comp].Fratio;
208 else
209 F = 1.0;
210
211 RochePot = Binary[Comp].RochePot;
212 RadiusBound = 0.8*Binary[Comp].Radius;
213 etacorr_over = Binary[Comp].LimEta;
214
215
216 /* ---------- Maximum Radius to search for Root ---------------- */
217
218 if (Flags.fill == OFF) {
219 RXCrit = BinPtr->RXCrit;
220 } else {
221 RXCrit = BinPtr->LimRad
222 + (0.07*SQR(Binary[Primary].Mq))
223 + 0.02*Binary[Primary].Mq;
224
225 if ( (Binary[Primary].Mq - 1.0) >= FLT_EPSILON)
226 RXCrit = BinPtr->LimRad
227 + (0.07*SQR(Binary[Secondary].Mq))
228 + 0.02*Binary[Secondary].Mq;
229
230 if ( (Binary[Primary].RocheFill - 1.05) <= FLT_EPSILON
231 && (Binary[Primary].Mq - 0.2) >= FLT_EPSILON
232 && (Binary[Primary].Mq - 5.0) <= FLT_EPSILON )
233 RXCrit = BinPtr->RXCrit + Binary[Primary].RocheFill - 1.;
234
235 else if ( (Binary[Primary].RocheFill - 1.02) <= FLT_EPSILON )
236 RXCrit = BinPtr->RXCrit + Binary[Primary].RocheFill - 1.;
237
238
239 /* critical angle w/r to z */
240 etacorr_over = BinPtr->RZatL1 /
241 sqrt( BinPtr->RZatL1 * BinPtr->RZatL1
242 + BinPtr->RLag1 * BinPtr->RLag1);
243 etacorr_over = asin( CLAMP(etacorr_over, -1.0, 1.0) );
244 }
245
246
247 BinPtr->Surface = 0.0;
248 BinPtr->Gravity = 0.0;
249
250 if (Comp == 0) {
251 T1 = BinPtr->Temperature;
252 T2 = Binary[Comp+1].Temperature;
253 } else {
254 T1 = BinPtr->Temperature;
255 T2 = Binary[Comp-1].Temperature;
256 }
257 Tr = T2/T1;
258 Tr = Tr * Tr * Tr * Tr;
259 n_eta = StepsPerPi;
260 n_phi = 0;
261
262 if (Flags.fill == OFF) BinPtr->LimEta = 0.0;
263
264 delta_eta = (PI - BinPtr->LimEta)/n_eta;
265 eta = PI-(0.5*delta_eta); /* start value, point opposite to companion */
266
267 if (Flags.fill == ON) { /* Throat adjustment */
268 etacorr_over = (BinPtr->LimEta - etacorr_over)
269 + 3.0 * delta_eta;
270 }
271
272 NumElem = 0;
273
274 /* ------------------- loop over eta ------------------------------ */
275
276 BinPtr->N_Rings = n_eta;
277
278 for (i = 0; i < n_eta; ++i) {
279
280 eta_corr = eta; /* Throat adjustment */
281 deltaeta_corr = delta_eta;
282
283 phi = 0.0;
284 sineta = sin(eta);
285 coseta = cos(eta);
286 lambda = coseta;
287
288
289 /* fix n_phi for GL animation, because the vertex list needs a */
290 /* fixed number vertices per phi loop */
291 #ifdef _WITH_OPENGL
292 if (Flags.animate == ON && Flags.use_opengl == ON) {
293 /* n_phi fixed to a power of 2 for GL animation */
294 n_phi = 32;
295 } else {
296 /* n_phi proportional to sin(eta) */
297 n_phi = (long) (10 + 1 + (int)(fabs(2 * n_eta * sineta)));
298 }
299 #else
300 /* n_phi proportional to sin(eta) */
301 n_phi = (long) (10 + 1 + (int)(fabs(2 * n_eta * sineta)));
302 #endif
303
304 delta_phi = (PI+PI)/n_phi;
305
306 /* check if enough memory */
307 if ((NumElem+n_phi) >= (int) MaximumElements) {
308 #ifdef _WITH_GTK
309 if (Flags.interactive == OFF)
310 nf_error(_(errmsg[10]));
311 else {
312 if (Flags.WantFit == OFF && Flags.WantMap == OFF)
313 make_dialog(_(errmsg[10]));
314 return(8);
315 }
316 #else
317 nf_error(_(errmsg[10]));
318 #endif
319 }
320
321 BinPtr->N_PhiStep[i] = n_phi; /* store number of phi intervals */
322
323 /* --------- Loop over phi ---------------------------------------- */
324
325 for (j = 0; j < n_phi; ++j) {
326
327 sinphi = sin(phi);
328 cosphi = cos(phi);
329
330 /* Throat adjustment for last three rings */
331 /* The throat is wider in x-y plane, thus we need to adjust eta, */
332 /* otherwise there would be an increasing gap towards higher */
333 /* latitudes. */
334
335 if ((Flags.fill == ON) && (i > (n_eta-4))) {
336 if (i == n_eta-3) eta_corr = eta
337 - ((etacorr_over/6.0)
338 -(delta_eta/2.0))*fabs(sinphi);
339 else if (i == n_eta-2) eta_corr = eta
340 - ((etacorr_over/2.0)
341 -(1.5*delta_eta))*fabs(sinphi);
342 else if (i == n_eta-1) eta_corr = eta
343 - ((5.0*etacorr_over/6.0)
344 -(5.0*delta_eta/2.0))*fabs(sinphi);
345
346 sineta = sin(eta_corr);
347 coseta = cos(eta_corr);
348 lambda = coseta;
349 deltaeta_corr = delta_eta+((etacorr_over/3.0)-delta_eta)*fabs(sinphi);
350 }
351
352 mu = sineta*cosphi;
353 nu = sineta*sinphi;
354
355 rad1 = RootFind(RocheSurface, RadiusBound, RXCrit,
356 100.0 * DBL_EPSILON, /* was 1.0e-8, */
357 "DivideSurface1", &testerr);
358 if (testerr == 1) return(8);
359
360 sqr_rad1 = rad1*rad1;
361 rho = 1.0/sqrt(1.0+(sqr_rad1)-2.0*rad1*lambda);
362 cub_rho = rho*rho*rho;
363
364 /* old code
365 C_rad1 = -1.0/(sqr_rad1) + Mq*(cub_rho*(lambda-rad1)-lambda)
366 +F*F*(Mq+1.0)*(sqr_rad1)*(1.0- nu*nu);
367 */
368 /* 2002-05-03 rwichmann error in expression fixed
369 */
370 C_rad1 = -1.0/(sqr_rad1) + Mq*(cub_rho*(lambda-rad1)-lambda)
371 +F*F*(Mq+1.0)*rad1*(1.0- nu*nu);
372 C_eta = sineta
373 *(Mq*(1.0-cub_rho)-(Mq+1.0)*rad1*F*F
374 *lambda*sinphi*sinphi);
375 C_phi = -(Mq + 1.0)*F*F*rad1*nu*cosphi;
376 C_x = C_rad1*lambda - C_eta*sineta;
377 C_y = C_rad1*mu + C_eta*lambda*cosphi - C_phi*sinphi;
378 C_z = C_rad1*nu + C_eta*lambda*sinphi + C_phi*cosphi;
379
380 grav = sqrt( C_x*C_x + C_y*C_y + C_z*C_z );
381
382 lvec = -C_x/grav;
383 mvec = -C_y/grav;
384 nvec = -C_z/grav;
385
386 DSurf = (sqr_rad1) * sineta * deltaeta_corr * delta_phi;
387
388 /* 2002-05-03 rwichmann HIGH_PRECISION code */
389
390 #ifdef HIGH_PRECISION
391 /*
392 * dA = sin(eta) deta dphi * ( R*R e(r) - R dR/deta e(eta)
393 * - R/sin(eta) dR/dphi e(phi) )
394 */
395
396 /* the dR/deta term
397 */
398 rdiff = NumDiffR (deltaeta_corr, eta_corr, sinphi,
399 RadiusBound, RXCrit, &testerr, eps);
400 if (testerr == 1) return (8);
401 rdiff = rad1 * rdiff * sineta * deltaeta_corr * delta_phi;
402
403 /* the dR/dphi term
404 */
405 rdiff2 = NumDiffRPhi (delta_phi, phi, sineta,
406 RadiusBound, RXCrit, &testerr, eps);
407 if (testerr == 1) return (8);
408 rdiff2 = rad1 * rdiff2 * deltaeta_corr * delta_phi;
409
410 /* absolute value
411 */
412 DSurf = sqrt(DSurf*DSurf + rdiff*rdiff + rdiff2*rdiff2);
413 #else
414 /*
415 * cos (e) = lambda*lvec+mu*mvec+nu*nvec
416 * = angle to surface normal
417 */
418 DSurf = DSurf/(lambda*lvec+mu*mvec+nu*nvec);
419 #endif
420
421 SurfPtr->eta = (float) eta_corr;
422 SurfPtr->phi = (float) phi;
423 SurfPtr->rad = (float) rad1;
424 SurfPtr->lambda = (float) (lambda * rad1);
425 SurfPtr->mu = (float) (mu * rad1);
426 SurfPtr->nu = (float) (nu * rad1);
427
428 SurfPtr->rho = (float) rho;
429 SurfPtr->grav = (float) grav;
430 SurfPtr->l = (float) lvec;
431 SurfPtr->m = (float) mvec;
432 SurfPtr->n = (float) nvec;
433
434
435 SurfPtr->area = (float) DSurf;
436
437 SurfPtr->MyNum = NumElem;
438 SurfPtr->ring = i;
439
440 Surf = Surf + DSurf;
441 MeanGrav = MeanGrav + (grav*DSurf);
442 MeanDarkGrav = MeanDarkGrav
443 + exp(BinPtr->GravDarkCoeff * log(grav))*DSurf;
444
445
446 /* --------------- INITIALIZE SURFACE TEMP ----------------- */
447
448 /* If we want to consider limb brightening for illuminated side
449 * of star, we need to mark illuminated/non-illuminated surface
450 * elements. Let the default be non-illuminated, and mark illuminated
451 * surface elements in LightReflect().
452 */
453 SurfPtr->AreaType = BACKSIDE;
454
455 SurfPtr->temp = T1;
456
457 phi = phi + delta_phi; /* increment phi */
458 ++NumElem; /* increment surface element count */
459 ++SurfPtr;
460
461 } /* ----------- end loop over phi --------------------------- */
462
463 eta = eta - delta_eta;
464
465 } /* -------------- end loop over eta ------------------------------- */
466
467
468 BinPtr->NumElem = NumElem;
469 BinPtr->NumPlus = NumElem;
470
471 if (Flags.debug[VERBOSE] == ON) {
472 printf("\n");
473 printf(_(" Divide Star %3d : %6ld Surface elements\n"),
474 Comp+1, BinPtr->NumElem);
475 }
476
477 BinPtr->Surface = Surf;
478 BinPtr->Gravity = MeanGrav/Surf;
479 BinPtr->DarkeningGrav = MeanDarkGrav/Surf;
480
481
482 /* ------------------- SET next/last_ptr ---------------------------- */
483
484 /* we need this only once */
485
486 if (Flags.first_pass == ON) {
487
488 SurfPtr = Surface[Comp];
489
490 n_last = 0;
491 n_next = 0;
492 NumElem = 0;
493
494 for (i = 0; i < n_eta; ++i) {
495
496 if (i > 0) n_last = n_phi; /* pointer setting */
497 n_phi = BinPtr->N_PhiStep[i];
498 if ( i < n_eta-1) n_next = BinPtr->N_PhiStep[i+1];
499
500 for (j = 0; j < n_phi; ++j) {
501
502
503 /* ----------- the first ring ---------------------------- */
504 if (i == 0) {
505 SurfPtr->last_ptr = NULL;
506 if (n_next == n_phi) {
507 SurfPtr->next_ptr =
508 &Surface[Comp][NumElem+n_next];
509 } else {
510 n_comp = ROUND(((float)(n_next)/(float)(n_phi))*(float)(j));
511 n_comp = n_phi - j + n_comp;
512 SurfPtr->next_ptr =
513 &Surface[Comp][NumElem+n_comp];
514 }
515 }
516 /* ----------- somewhere in the middle -------------------- */
517 else if ( (i > 0) && (i < (n_eta-1))) {
518 if (n_next == n_phi) {
519 NextPtr = NumElem+n_next;
520 if (NextPtr >= BinPtr->NumElem)
521 NextPtr = BinPtr->NumElem-1;
522 SurfPtr->next_ptr =
523 &Surface[Comp][NextPtr];
524 } else {
525 n_comp = ROUND(((float)(n_next)/(float)(n_phi))*(float)(j));
526 n_comp = n_phi - j + n_comp;
527 NextPtr = NumElem+n_comp;
528 if (NextPtr >= BinPtr->NumElem)
529 NextPtr = BinPtr->NumElem-1;
530 SurfPtr->next_ptr =
531 &Surface[Comp][NextPtr];
532 }
533 if (n_last == n_phi) {
534 SurfPtr->last_ptr =
535 &Surface[Comp][NumElem-n_next];
536 } else {
537 n_comp = ROUND(((float)(n_last)/(float)(n_phi))*(float)(j));
538 n_comp = (n_last - n_comp) + j;
539 SurfPtr->last_ptr =
540 &Surface[Comp][NumElem-n_comp];
541 }
542 }
543
544 /* ----------- the last ring ---------------------------- */
545
546 /*@ignore@*/
547 /* Storage SurfPtr->last_ptr is kept in one path, but live in */
548 /* another. */
549
550 else if (i == (n_eta-1)) {
551
552 if (n_last == n_phi) {
553 NextPtr = NumElem-n_phi;
554 if (NextPtr >= BinPtr->NumElem)
555 NextPtr = BinPtr->NumElem-1;
556 SurfPtr->last_ptr =
557 &Surface[Comp][NextPtr];
558 } else {
559 n_comp = ROUND(((float)(n_last)/(float)(n_phi))*(float)(j));
560 n_comp = (n_last - n_comp) + j;
561 NextPtr = NumElem-n_comp -1;
562 if (NextPtr >= BinPtr->NumElem)
563 NextPtr = BinPtr->NumElem-1;
564 SurfPtr->last_ptr =
565 &Surface[Comp][NextPtr];
566 }
567
568 if (Flags.fill == OFF) {
569 n_comp = (j + (n_phi/2)) % n_phi;
570 n_comp = - j + n_comp;
571 NextPtr = NumElem+n_comp;
572 SurfPtr->next_ptr =
573 &Surface[Comp][NextPtr];
574 } else {
575 if (Comp == Primary) {
576 SurfPtr->next_ptr =
577 &Surface[Secondary][NextPtr];
578 } else {
579 NextPtr = NumElem;
580 SurfPtr->next_ptr =
581 &Surface[Primary][NextPtr];
582 }
583 }
584 }
585
586
587 if (NumElem != 0) {
588 if (j != 0) {
589 SurfPtr->phi_ptr =
590 &Surface[Comp][NumElem-1];
591 } else {
592 SurfPtr->phi_ptr =
593 &Surface[Comp][NumElem+n_phi-1];
594 }
595 if (j != (n_phi-1)) {
596 SurfPtr->phi1_ptr =
597 &Surface[Comp][NumElem+1];
598 } else {
599 SurfPtr->phi1_ptr =
600 &Surface[Comp][NumElem-(n_phi-1)];
601 }
602 } else {
603 SurfPtr->phi1_ptr =
604 &Surface[Comp][1];
605 SurfPtr->phi_ptr =
606 &Surface[Comp][1];
607 }
608
609 ++NumElem;
610 ++SurfPtr;
611
612 }
613 /*@end@*/
614 }
615
616 }
617
618 /* -------------- END SET next/last_ptr ----------------------------- */
619
620 tcalc:
621
622 /**********************************************************************/
623 /* */
624 /* Temperature Distribution */
625 /* */
626 /**********************************************************************/
627
628 /* --------------- Gravitational Darkening ------------------------- */
629
630 /* T propto g**GravDarkCoeff */
631
632 SurfPtr = Surface[Comp];
633
634 for (i = 0; i < BinPtr->NumElem; ++i) {
635
636 SurfPtr->temp =
637 BinPtr->Temperature
638 * (exp(BinPtr->GravDarkCoeff * log(SurfPtr->grav))
639 / BinPtr->DarkeningGrav );
640
641 SurfPtr->orig_temp = SurfPtr->temp;
642
643 ++SurfPtr;
644
645 }
646
647 /* ------------------------ Spots -------------------- */
648
649
650 /* make spots or initialize the dimfactor to 1.0 */
651
652 if (Flags.first_pass == ON) {
653 if (Comp == Primary) {
654 if (Flags.Spots1 > OFF) {
655 MakeSpots(Primary, 0);
656 } else {
657 SurfPtr = Surface[Primary];
658 for (j = 0; j < Binary[Primary].NumElem; ++j) {
659 SurfPtr->SumDimFactor = 1.0;
660 ++SurfPtr;
661 }
662 }
663 }
664
665 if (Comp == Secondary) {
666 if (Flags.Spots2 > OFF) {
667 MakeSpots(Secondary, 0);
668 } else {
669 SurfPtr = Surface[Secondary];
670 for (j = 0; j < Binary[Secondary].NumElem; ++j) {
671 SurfPtr->SumDimFactor = 1.0;
672 ++SurfPtr;
673 }
674 }
675 }
676
677 } else {
678
679 /* move spots if asynchroneous rotation or elliptic orbit */
680
681 if (Comp == Primary) {
682 if ( (Flags.Spots1 > OFF) &&
683 ( Flags.asynchron1 == ON ||
684 Flags.asynchron2 == ON ||
685 Flags.elliptic == ON) )
686 MakeSpots(Primary, Orbit.PhaseIndex + N_FluxOutFull);
687 }
688
689 if (Comp == Secondary) {
690 if ( (Flags.Spots2 > OFF) &&
691 ( Flags.asynchron1 == ON ||
692 Flags.asynchron2 == ON ||
693 Flags.elliptic == ON) )
694 MakeSpots(Secondary, Orbit.PhaseIndex + N_FluxOutFull);
695 }
696
697 }
698
699 /* --------------------------------- Reflection -------------------- */
700
701 if ((Flags.reflect > OFF) && (Comp == Secondary)) {
702 LightReflect();
703 }
704
705 if (Flags.reflect == OFF) SimpleReflect(Comp);
706
707 /* ---------------------------- Normalize Power -------------------- */
708 /* (for testing only) */
709
710 #if 0
711 if (Flags.elliptic == ON)
712 {
713 double Power = 0.0;
714 double CFactor;
715 double CPower;
716 double pt;
717 SurfPtr = Surface[Comp];
718
719 for (i = 0; i < BinPtr->NumElem; ++i)
720 {
721 pt = SurfPtr->temp;
722 Power += (pt * pt) * (pt * pt) * SurfPtr->area;
723 ++SurfPtr;
724 }
725
726 if (NormPower[Comp] < FLT_EPSILON)
727 {
728 NormPower[Comp] = Power;
729
730 fprintf(stderr, "#FIXME <%8.5e> 0.0 %d\n", Power, Comp);
731 }
732 else
733 {
734 CFactor = NormPower[Comp]/Power;
735 CFactor = sqrt(sqrt(CFactor));
736
737 SurfPtr = Surface[Comp];
738 CPower = 0.0;
739
740 for (i = 0; i < BinPtr->NumElem; ++i)
741 {
742 SurfPtr->temp *= CFactor;
743 pt = SurfPtr->temp;
744 CPower += (pt * pt) * (pt * pt) * SurfPtr->area;
745 ++SurfPtr;
746 }
747 fprintf(stderr, "#FIXME <%8.5e> <%8.5e> %8.5e\n",
748 Power, CPower, CFactor);
749 }
750 }
751 #endif
752
753 /* ---------------------------- mean temperature ------------------- */
754
755 Weight = 0.0; TempMean = 0.0;
756 SurfPtr = Surface[Comp];
757 for (i = 0; i < BinPtr->NumElem; ++i) {
758 Area = SurfPtr->area;
759 TempMean = TempMean + SurfPtr->temp * Area;
760 Weight = Weight + Area;
761 /*------- Sollte hier nicht irgendwo ein SurfPtr++ stehen????? */
762 /* eigentlich schon ... RW Tue Sep 11 12:54:20 CEST 2001 */
763 ++SurfPtr;
764 }
765 BinPtr->TempMean = TempMean/Weight;
766
767
768
769 /**********************************************************************/
770 /* */
771 /* Flux intensity computation */
772 /* */
773 /**********************************************************************/
774
775 if (Comp == Secondary) {
776 if (Flags.blackbody == ON) {
777 BlackbodyFlux(Primary);
778 BlackbodyFlux(Secondary);
779 } else {
780 ModelFlux(Primary);
781 ModelFlux(Secondary);
782 }
783 }
784
785 /* ---------------- output --------------------- */
786
787 if (Flags.debug[VERBOSE] == ON) {
788 /* printf("\n Star: %5d\n", (Comp+1)); */
789 printf(_(" Surface Area, Mean Gravity (dimensionless): %8.6f %8.6f\n"),
790 BinPtr->Surface, BinPtr->Gravity);
791 printf(_(" Temperature is %8.6f\n"), BinPtr->Temperature);
792
793 printf("\n");
794 printf(_(" Temperature Distribution corrected for: \n"));
795 printf(_(" Reflection and Gravity Darkening\n"));
796 printf(_(" Coefficient for Gravity Darkening was: %6.4f\n"),
797 BinPtr->GravDarkCoeff);
798 }
799
800 /* ---------------- statistics --------------------- */
801
802 if (Comp == Secondary && Flags.debug[SURFACE] == ON)
803 {
804 for (i = 0; i < 2; ++i)
805 {
806 int ring;
807 int NumElem = Binary[i].NumElem;
808 double r_temp[STEPS_PER_PI] = { 0.0 };
809 double r_flux[STEPS_PER_PI] = { 0.0 };
810 double r_refl[STEPS_PER_PI] = { 0.0 };
811 int r_n[STEPS_PER_PI] = { 0 };
812
813 for (j = 0; j < NumElem; ++j)
814 {
815 ring = Surface[i][j].ring;
816 r_temp[ring] += Surface[i][j].orig_temp;
817 r_refl[ring] += Surface[i][j].temp;
818 r_flux[ring] += Surface[i][j].f_[Rmag];
819 ++r_n[ring];
820 }
821
822 printf("\n");
823 printf("%02ld Temp Trefl Ratio Flux\n",
824 i+1);
825 printf("------------------------------------------------------\n");
826 for (j = 0; j < Binary[i].N_Rings; ++j)
827 {
828 printf("%02ld %10.4g %10.4g %10.4g %10.4g %03d\n", j,
829 r_temp[j]/r_n[j],
830 r_refl[j]/r_n[j],
831 r_refl[j]/r_temp[j],
832 r_flux[j]/r_n[j],
833 r_n[j]);
834 }
835 printf("\n");
836 }
837 }
838 return(0);
839
840 }
841
842 #if 0
843 /* Can we rename this to a more intuitive name light rotatex or what ever;
844 I don't know what "tile" in this context means; Why using floats and
845 not double values ? MK */
846
847 /* This is a no-op if the tilt angle is zero
848 */
849 static void TileDisk1(double TiltAngle,SurfaceElement *SurfPtr)
850 {
851 float lambda_strich;
852 float nu_strich;
853
854 lambda_strich = (float) SurfPtr -> lambda*cos(TiltAngle)-
855 SurfPtr -> nu*sin(TiltAngle);
856 nu_strich = (float) SurfPtr -> lambda*sin(TiltAngle)+
857 SurfPtr -> nu*cos(TiltAngle);
858
859 SurfPtr->lambda = lambda_strich;
860 SurfPtr->nu = nu_strich;
861 }
862
863 /* This is a no-op if the tilt angle is zero
864 */
865 static void TileDisk2(double TiltAngle,SurfaceElement *SurfPtr)
866 {
867 float l_strich;
868 float n_strich;
869
870 l_strich = (float) SurfPtr -> l*cos(TiltAngle)+
871 SurfPtr -> n*sin(TiltAngle);
872 n_strich = (float) SurfPtr -> l*sin(TiltAngle)-
873 SurfPtr -> n*cos(TiltAngle);
874
875 SurfPtr -> l = l_strich;
876 SurfPtr -> n = n_strich;
877 }
878
879 /* This is a no-op if the tilt angle is zero
880 */
881 static void TileDisk3(double TiltAngle,SurfaceElement *SurfPtr)
882 {
883 float lambda_strich;
884 float mu_strich;
885
886 lambda_strich = (float) SurfPtr -> lambda*cos(TiltAngle)-
887 SurfPtr -> mu*sin(TiltAngle);
888
889 mu_strich = (float) SurfPtr -> lambda*sin(TiltAngle)-
890 SurfPtr -> mu*cos(TiltAngle);
891
892 SurfPtr -> lambda = lambda_strich;
893 SurfPtr -> mu = mu_strich;
894 }
895
896 static void CalculateLayerVecE1(SurfaceElement *SurfPtr)
897 {
898 double vec_en[3]; /* Normalenvektor der Ebene */
899
900 vec_en[0] = SurfPtr->l;
901 vec_en[1] = SurfPtr->m;
902 vec_en[2] = SurfPtr->n;
903
904 if (vec_en[0] != 0)
905 {
906 SurfPtr->VecE1[0] = -((vec_en[1]+vec_en[2])/vec_en[0]);
907 SurfPtr->VecE1[1] = 1.0;
908 SurfPtr->VecE1[2] = 1.0;
909 }
910 else if (vec_en[1] != 0)
911 {
912 SurfPtr->VecE1[0] = 1.0;
913 SurfPtr->VecE1[1] = -((vec_en[0]+vec_en[2])/vec_en[1]);
914 SurfPtr->VecE1[2] = 1.0;
915 }
916 else
917 {
918 SurfPtr->VecE1[0] = 1.0;
919 SurfPtr->VecE1[1] = 1.0;
920 SurfPtr->VecE1[2] = -((vec_en[0]+vec_en[1])/vec_en[2]);
921 }
922 }
923
924 static void CalculateLayerVecE2(SurfaceElement *SurfPtr)
925 {
926 double vec_en[3]; /* Normalenvektor der Ebene */
927 double vec_e1[3]; /* Vektor Ebene 1 */
928
929 vec_en[0] = SurfPtr->l;
930 vec_en[1] = SurfPtr->m;
931 vec_en[2] = SurfPtr->n;
932
933 vec_e1[0] = SurfPtr->VecE1[0];
934 vec_e1[1] = SurfPtr->VecE1[1];
935 vec_e1[2] = SurfPtr->VecE1[2];
936
937 SurfPtr->VecE2[0] = vec_en[1]*vec_e1[2]-vec_en[2]*vec_e1[1];
938 SurfPtr->VecE2[1] = vec_en[2]*vec_e1[0]-vec_en[0]*vec_e1[2];
939 SurfPtr->VecE2[2] = vec_en[0]*vec_e1[1]-vec_en[1]*vec_e1[0];
940
941 }
942 #endif
943
944 #ifdef HAVE_DISK
945
946 /* This macro yields correct results only if a > 0 and b in [0, 2*PI).
947 */
948 #define DIST_CIRCLE(a, b) fabs(PI - fabs(PI - fmod((a) - (b), (2.0 * PI) ) ))
949
950
951 /****************************************************************************
952 @author Rainer Wichmann
953 @short Compute height of hot spot
954 @long Compute the height of the hot spot depending on distance to centre.
955 @param (double) Maximum height (at zero distance)
956 @param (double) Normalized distance (0..1) to hot spot centre
957 @return (double) Height at given distance
958 @heading Disk Setup
959 ****************************************************************************/
puffHeight(double h0,double dist)960 double puffHeight(double h0, double dist)
961 {
962 double ret = h0 * /* exp(-(dist)) */ cos (dist * (PI/2.0));
963 return ret;
964 }
965
966 /****************************************************************************
967 @author Rainer Wichmann
968 @short Computed puffed hot spot
969 @long Computed puffed hot spot
970 @param (SurfaceElement*) An individual surface element in the spot area
971 @param (double) Normalized distance (0..1) to hot spot centre
972 @param (double) Outer radius
973 @param (double) Position along circumference
974 @param (double) inner height before puffing
975 @param (double) outer height before puffing
976 @param (double) Surface area not adjusted for cos(rise_angle)
977 @return (double) Adjusted surface area of element
978 @heading Disk Setup
979 ****************************************************************************/
puff_element(SurfaceElement * SurfPtr,double circ_dist,double Rout,double phi,double h_i_old,double h_i_new,double DSurf)980 double puff_element(SurfaceElement *SurfPtr,
981 double circ_dist,
982 double Rout,
983 double phi,
984 double h_i_old,
985 double h_i_new,
986 double DSurf /* to compute corrected area */)
987 {
988 /* need to re-compute:
989 * height SurfPtr->nu,
990 * normal vector SurfPtr->l, SurfPtr->m, SurfPtr->n
991 * area SurfPtr->area
992 */
993 double hs_puff = Binary[Disk].puffHS; /* puffyness */
994 double hs_depth = Binary[Disk].depthHS; /* radial depth */
995 double DSurfCor = DSurf; /* surface of element, sin included */
996
997 double RLag1 = Binary[Secondary].RLag1;
998
999 double radial_dist;
1000 double puff = 0.0;
1001
1002 hs_puff *= RLag1;
1003 hs_depth *= RLag1;
1004
1005 if (SurfPtr->AreaType == BOTTOM_SEGMENT || SurfPtr->AreaType == TOP_SEGMENT)
1006 {
1007 double height_inner = h_i_old;
1008 double height_outer = h_i_new;
1009 double rise_angle = 1.0;
1010 double sign = 1.0;
1011
1012 double radial_dist_i;
1013 double radial_dist_o;
1014 double puff_o = 0.0;
1015 double puff_i = 0.0;
1016
1017 /* avoid division by zero
1018 */
1019 if (hs_depth > DBL_EPSILON)
1020 {
1021 hs_puff = puffHeight(hs_puff, circ_dist);
1022
1023 /* No need to safeguard here, because this function is only called
1024 * for disk regions within the hot spot.
1025 */
1026 radial_dist = (SurfPtr->rad - (Rout - hs_depth))/hs_depth;
1027 radial_dist_i = (SurfPtr->RadiusIn - (Rout - hs_depth))/hs_depth;
1028 radial_dist_o = (SurfPtr->RadiusOut - (Rout - hs_depth))/hs_depth;
1029
1030 puff = hs_puff * radial_dist * radial_dist;
1031 puff_o = hs_puff * radial_dist_o * radial_dist_o;
1032 puff_i = (radial_dist_i > 0.0) ? (hs_puff * radial_dist_i * radial_dist_i) : 0.0;
1033 }
1034
1035 /* invert sign for bottom part of disk
1036 */
1037 if (SurfPtr->AreaType == BOTTOM_SEGMENT)
1038 {
1039 puff *= (-1.0);
1040 puff_o *= (-1.0);
1041 puff_i *= (-1.0);
1042 rise_angle = (-1.0);
1043 sign = (-1.0);
1044 }
1045
1046 /* Corrected height
1047 */
1048 SurfPtr->nu += puff;
1049
1050 /* Corrected surface area
1051 */
1052 height_inner += puff_i;
1053 height_outer += puff_o;
1054
1055 /* This is ok because rise_angle is initialized to +/- 1.0
1056 */
1057 rise_angle *= atan( fabs(height_outer - height_inner) / fabs(SurfPtr->RadiusOut - SurfPtr->RadiusIn) );
1058 DSurfCor /= cos( rise_angle );
1059 SurfPtr->area = DSurfCor;
1060
1061 /* Corrected normal vector
1062 */
1063 SurfPtr->l = (float) (-cos(phi)) * sin(rise_angle) * sign;
1064 if (fabs(SurfPtr->l) <= DBL_EPSILON) {
1065 SurfPtr->l = 0.0;
1066 }
1067 SurfPtr->m = (float) (-sin(phi)) * sin(rise_angle) * sign;
1068 if (fabs(SurfPtr->m) <= DBL_EPSILON) {
1069 SurfPtr->m = 0.0;
1070 }
1071 SurfPtr->n = (float) cos(rise_angle) * sign;
1072 if (fabs(SurfPtr->n) <= DBL_EPSILON) {
1073 SurfPtr->n = 0.0;
1074 }
1075
1076 }
1077 else if (SurfPtr->AreaType == OUT_RECTANGLE)
1078 {
1079 double height_max = h_i_old;
1080 double height_this = h_i_new;
1081 double sfactor = 1.0;
1082
1083 /* avoid division by zero
1084 */
1085 if (hs_depth > DBL_EPSILON && height_max > DBL_EPSILON)
1086 {
1087 hs_puff = puffHeight(hs_puff, circ_dist);
1088
1089 sfactor = hs_puff / height_max;
1090 puff = sfactor * height_this;
1091
1092 DSurfCor *= sfactor;
1093 SurfPtr->AreaHeight *= sfactor;
1094 SurfPtr->nu += puff; /* Fixed 2010-02-10, = -> += */
1095 }
1096 }
1097
1098 /* Return area for addition to total area
1099 */
1100 return DSurfCor;
1101 }
1102
1103 /****************************************************************************
1104 @author Patrick Risse
1105 @short Divide disk surface into grid
1106 @long First, set up the grid, then compute the
1107 temperature distribution on the grid, including
1108 spots and reflection
1109 @param (int) Comp The stellar component
1110 @return (int) The exit status
1111 @heading Disk Setup
1112 ****************************************************************************/
DivideDisk(int Comp)1113 int DivideDisk(int Comp)
1114 {
1115 register long i,j; /* loop variables */
1116 int n_phi = 32.0; /* surface steps in phi */
1117 long NumElem; /* current surface element */
1118 /*long NextPtr;*/ /* next surface element */
1119 /*int testerr = 0;*/ /* exit status of RootFind */
1120
1121 double phi; /* phi of surface element */
1122 double T1; /* temperature of this star */
1123 double Rout,Rin; /* Outer/inner Radius of the Disk */
1124 #if 0
1125 double RadiusBound = 0.0; /* inner bound for RootFind */
1126 double R1; /* Roche fill of this star */
1127 double etacorr_over; /* Throat adjustment */
1128 #endif
1129 /*long n_next, n_last, n_comp;*/ /* Pointer setting */
1130 /* for mean temperature calculation */
1131 double Area; /* area of element */
1132
1133 double Weight=0.0; /* total area */
1134 double TempMean=0.0; /* mean temperature */
1135 double Phase; /* Phase of the System. Is needed
1136 for recalculating the Disk Position*/
1137
1138 SurfaceElement *SurfPtr; /* pointer to surface */
1139 BinaryComponent *BinPtr; /* pointer to binary */
1140
1141 /* ---------------- Variables added for the disk by P.R. ---------------*/
1142 double Thick_out; /* Thickness of the disk on the
1143 outside */
1144 double Thick_in; /* Thickness of the disk on the
1145 inside */
1146 double diskExp = 1.0; /* Exponent for r dependency of
1147 scale height */
1148 double diskProp = 1.0; /* Proportionality factor for disk
1149 height */
1150 double diskAdd = 1.0; /* Additive constant for disk
1151 height */
1152 double d_s1d; /* Area of one element on the top/
1153 bottom side of the disk */
1154 double d_s2d_out; /* Area */
1155 double d_s2d_in;
1156 double P_0;
1157 double DSurf; /* surface of element */
1158 double DSurfCor; /* surface of element, sin included */
1159 double Surf = 0.0; /* total surface */
1160 double ring_height;
1161 double r_i_old;
1162 double r_i_new;
1163 double h_i_old;
1164 double h_i_new;
1165 double rise_angle;
1166 double r_inner_sqr;
1167 double r_outer_sqr;
1168 double dphi = 0.0;
1169 double r_j=0;
1170 double alpha; /*Angles in Degree*/
1171
1172 double offset_in=0;
1173 double offset_out=0;
1174 int n_r; /* Number of Rings on the top/bottom
1175 of the disk */
1176 /*long TopElem,OutElem,BottomElem,InElem;*/
1177
1178 int NumOfTopElements,NumOfOuterElements ;
1179 int NumOfBottomElements,NumOfInnerElements;
1180 int NumOfTotalElements;
1181
1182 double warpfactor[42]; /* Inclination of each ring for the warp*/
1183 #if 0
1184 float rot[42];
1185 int warp_no; /* Entry_No of warpfactor */
1186 int tilt_flag = 0;
1187 double TiltAngle;
1188 #endif
1189
1190 double t_c;
1191 double t_temp;
1192 int t_flag = Flags.tdisk;
1193
1194 double sin_phi, cos_phi, sin_alpha, cos_alpha;
1195
1196 double hs_long = internal_longitude(); /* spot longitude */
1197 double hs_ext = ((Binary[Disk].extentHS / 2.0) * DTOR); /* extent in phi */
1198 double hs_depth = Binary[Disk].depthHS; /* depth in R */
1199 double hs_temp = Binary[Disk].tempHS; /* temperature */
1200 double hs_puff = Binary[Disk].puffHS; /* puffyness */
1201 double hs_phi; /* temp variable */
1202 double hs_dist; /* temp variable */
1203
1204 Binary[Disk].longitude_intern = hs_long;
1205
1206 n_r = floor(StepsPerPi/1.5);
1207
1208 if (t_flag == 0) /* simple */
1209 diskExp = 1.0;
1210 else if (t_flag == 1) /* isothermal */
1211 diskExp = 1.5;
1212 else if (t_flag == 2) /* reprocessing */
1213 diskExp = (9.0/8.0);
1214
1215 /*BinPtr->N_PhiStep[i] = n_phi;*/ /* store number of phi intervals */
1216
1217 /* ---------- initialize ---------------------------------------------- */
1218
1219 SurfPtr = Surface[Comp];
1220 BinPtr = &Binary[Comp];
1221
1222 Mq = Binary[Comp].Mq;
1223 if (Comp == Primary && Flags.asynchron1 == ON)
1224 F = Binary[Comp].Fratio;
1225 else if (Comp == Secondary && Flags.asynchron2 == ON)
1226 F = Binary[Comp].Fratio;
1227 else
1228 F = 1.0;
1229
1230 /* RLag1 = Binary[Comp].RLag1;
1231 * This Value is not used. Why is it defined? P.R. 19.6.2001
1232 */
1233 RochePot = Binary[Comp].RochePot;
1234 /* RadiusBound = 0.8*Binary[Comp].Radius; */
1235 /* etacorr_over = Binary[Comp].LimEta; */
1236
1237 Phase = 0.0; /* ((2*PI/PhaseSteps) * index) */;
1238
1239 BinPtr->RocheFill = Binary[Secondary].RocheFill;
1240 BinPtr->RCrit = Binary[Secondary].RCrit;
1241 BinPtr->RLag1 = Binary[Secondary].RLag1;
1242
1243 /* RLag1 is not calculated for the Disk so i take it from the
1244 * Secondary Star. There must be a better possibility to define
1245 * the RLag1 for the Disk. ToDo: Where is RLag1 defined?
1246 * P.R. 19.06.2001
1247 */
1248
1249 Rout = BinPtr->RLag1 * BinPtr->Rout;
1250 Rin = BinPtr->RLag1 * BinPtr->Rin; /*BinPtr->RocheFill; */
1251
1252 hs_puff = hs_puff * BinPtr->RLag1;
1253 hs_depth = hs_depth * BinPtr->RLag1;
1254
1255 BinPtr->Surface = 0.0;
1256 BinPtr->Gravity = 0.0;
1257
1258 T1 = BinPtr->Temperature;
1259 /* R1 = BinPtr->Radius; */
1260
1261 NumElem = 0;
1262
1263 BinPtr->Segments = n_phi;
1264
1265 /* We need to define the proportionality factor such that the */
1266 /* disk temperature matches at Rin: T1 = t_c * pow(Rin, -0.75) */
1267
1268 t_c = (T1 / pow(Rin, -0.75));
1269
1270 /* We need to define the proportionality factor such that the disk */
1271 /* height matches at Rin: Thick = diskProp * pow(Rin, diskExp) */
1272
1273 if (t_flag == 0)
1274 {
1275 diskAdd = BinPtr->Thick;
1276 diskProp = BinPtr->HR;
1277 }
1278 else
1279 {
1280 diskAdd = 0.0;
1281 diskProp = (BinPtr->Thick / pow(Rin, diskExp));
1282 }
1283
1284 /* damn is this a hard job ! It's difficult to set the number of */
1285 /* phi steps to n_phi, this finally took 1 year to be done MK */
1286 /* BinPtr->N_PhiStep[0] = n_phi; */
1287
1288 /* This Function is not tested. There are some bugs inside */
1289 /* I will fix it later MPR 12.09.2001 */
1290
1291 if (Flags.warp == ON ) {
1292 /* Create warpfactor */
1293 printf("Creat standard warp\n");
1294 for (i=0; i< n_r-1; i++)
1295 {
1296 if (i <= 15) warpfactor[i]=(15.0/360)*2*PI;
1297 if (i > 15) warpfactor[i]=((15.0 -((i-15)*0.75))/360)*2*PI;
1298 if (warpfactor[i] <= 0) warpfactor[i]=0.0;
1299 /*
1300 if (i >= 1)
1301 rot[i]=((2.5/360)*2*PI)*i;
1302 */
1303 }
1304 }
1305
1306
1307 /* ------------------- Start building Disk grid --------------------- */
1308
1309 /* --------------------- General Calculations ------------------------ */
1310
1311 /* HR := Height/Radius
1312 * phi := angle in xy-plane
1313 * n_phi := no. of sectors
1314 * n_r := no. of rings
1315 */
1316
1317 /* Thickness of the Disk on the outer side
1318 */
1319 Thick_out = diskAdd + diskProp * pow(Rout, diskExp);
1320
1321 /* Thickness of the Disk on the inner side
1322 */
1323 Thick_in = diskAdd + diskProp * pow(Rin, diskExp);
1324
1325 /* Average area of one element
1326 */
1327 d_s1d = (PI*((Rout*Rout)-(Rin*Rin)))/(n_r*n_phi);
1328
1329 /* Area of one element on the outer side
1330 */
1331 d_s2d_out = (2.0*PI*Rout*Thick_out)/n_phi;
1332
1333 /* Area of one element on the inner side
1334 */
1335 d_s2d_in = (2.0*PI*Rin*Thick_in)/n_phi;
1336
1337 /* Number of rings rounded (to make elements approx. equal sized)
1338 */
1339 BinPtr->RingsOut = floor((d_s2d_out/d_s1d)+0.5);
1340
1341 /* Number of rings rounded (to make elements approx. equal sized)
1342 */
1343 BinPtr->RingsIn = floor((d_s2d_in/d_s1d)+0.5);
1344
1345 P_0 = (Rout*Rout-Rin*Rin)/n_r;
1346 r_i_old = Rin;
1347
1348 /*
1349 TiltAngle = (BinPtr->Tilt/360)*2*PI;
1350 if (fabs(TiltAngle) > FLT_EPSILON)
1351 tilt_flag = 1;
1352 */
1353
1354 if (Flags.debug[GEO]) {
1355 printf("\nDisk Geometry\n----------------------------------------\n");
1356 printf("Sectors : %d \n", n_phi);
1357 printf("Rings : %d \n", n_r);
1358 printf("RingsOut: %d \n", BinPtr->RingsOut);
1359 printf("ThickOut: %f \n", Thick_out);
1360 printf("RingsIn : %d \n", BinPtr->RingsIn);
1361 printf("ThickIn : %f \n", Thick_in);
1362 }
1363
1364 /* ----------------- Calculate number of elements for the Disk------ */
1365
1366 NumOfTopElements = n_phi * n_r;
1367
1368 NumOfOuterElements = BinPtr->RingsOut * n_phi;
1369
1370 NumOfBottomElements = n_phi * n_r;
1371
1372 NumOfInnerElements = BinPtr->RingsIn * n_phi;
1373
1374 NumOfTotalElements = NumOfTopElements + NumOfOuterElements +
1375 NumOfBottomElements + NumOfInnerElements;
1376
1377 if (Flags.debug[GEO]) {
1378 printf("Elements: %d\n", NumOfTotalElements);
1379 }
1380
1381 /**********************************************************************/
1382 /* */
1383 /* ---------------------- Calculate top of disk -------------------- */
1384 /* */
1385 /**********************************************************************/
1386
1387 r_inner_sqr = (r_i_old*r_i_old);
1388
1389 /* Area of an element. This can be lifted outside the loop, since
1390 * the progression of radii ensures equal areas for each element.
1391 */
1392 DSurf = (PI*P_0)/n_phi;
1393
1394 /* ----------------- loop over rings (n_r) ---------------------- */
1395
1396 for (i = 0; i < n_r; i++) {
1397
1398 phi = 0.0;
1399
1400 /* check if enough memory */
1401 if ((NumElem+n_phi) >= (int) MaximumElements) {
1402
1403 if (Flags.debug[GEO])
1404 printf("Top of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements);
1405
1406 #ifdef _WITH_GTK
1407 if (Flags.interactive == OFF)
1408 nf_error(_(errmsg[10]));
1409 else {
1410 if (Flags.WantFit == OFF && Flags.WantMap == OFF)
1411 make_dialog(_(errmsg[10]));
1412 return(8);
1413 }
1414 #else
1415 nf_error(_(errmsg[10]));
1416 #endif
1417 }
1418
1419 /* Radius of next circle. This progression ensures that all rings
1420 * have equal surface if the rings are flat.
1421 */
1422 r_i_new = sqrt(r_inner_sqr + ((i+1)*P_0));
1423
1424 /* Center of element
1425 */
1426 r_j = (r_i_new+r_i_old)/2.0;
1427
1428 /* Height of this ring
1429 */
1430 ring_height = ( (diskAdd + diskProp * pow(r_j, diskExp)) / 2.0);
1431
1432 /* angle of rise for this ring
1433 */
1434 h_i_old = ( (diskAdd + diskProp * pow(r_i_old, diskExp)) / 2.0);
1435 h_i_new = ( (diskAdd + diskProp * pow(r_i_new, diskExp)) / 2.0);
1436
1437 rise_angle = atan ( fabs(h_i_new - h_i_old) / fabs(r_i_new - r_i_old) );
1438
1439 /* Corrected surface area
1440 */
1441 DSurfCor = DSurf / cos ( rise_angle );
1442
1443 /* shift rings in phi by random amount to reduce numerical artifacts
1444 */
1445 #ifdef _WITH_OPENGL
1446 if (Flags.use_opengl == OFF)
1447 #endif
1448 dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0)));
1449
1450 /* Angle of rise from the middle to the edge of the Accretion Disk
1451 */
1452 if (t_flag == 0)
1453 alpha = 0.5 * diskProp;
1454 else
1455 alpha = rise_angle /* 0.5 * diskProp * diskExp * pow(r_j, (diskExp - 1.0)) */;
1456
1457 sin_alpha = sin(alpha);
1458 cos_alpha = cos(alpha);
1459
1460 if (t_flag < 2) /* isothermal */
1461 t_temp = T1;
1462 else
1463 t_temp = t_c * pow(r_j, -0.75);
1464
1465 if (Flags.debug[GEO]) {
1466 printf("Ring : %ld \n", i);
1467 printf("alpha : %f \n", alpha);
1468 printf("height : %f \n", ring_height);
1469 printf("temp : %f \n", t_temp);
1470 }
1471
1472 /* ----------------- loop over phi -------------------------------- */
1473
1474 for (j = 0; j < n_phi; ++j) {
1475
1476 phi = ((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */
1477
1478 phi += dphi;
1479
1480 sin_phi = sin(phi);
1481 cos_phi = cos(phi);
1482
1483 /* PR added 07.05.02 is needed for eclipsing
1484 * by disk verification
1485 */
1486 SurfPtr->RadiusIn = r_i_old;
1487 SurfPtr->RadiusOut = r_i_new;
1488
1489 /* phi counted from x to y axis
1490 */
1491 SurfPtr->lambda = (float) r_j * cos_phi;
1492 SurfPtr->mu = (float) r_j * sin_phi;
1493
1494 /* z is thickness/2.0 + (H/R * r_j)/2.0
1495 */
1496 SurfPtr->nu = (float) ring_height;
1497
1498 /* for a flaring disk, the normal vector points up and inwards,
1499 * so we need to switch sign for l, m
1500 */
1501 SurfPtr->l = (float) (-cos_phi) * sin_alpha;
1502 if (fabs(SurfPtr->l) <= DBL_EPSILON) {
1503 SurfPtr->l = 0.0;
1504 }
1505 SurfPtr->m = (float) (-sin_phi) * sin_alpha;
1506 if (fabs(SurfPtr->m) <= DBL_EPSILON) {
1507 SurfPtr->m = 0.0;
1508 }
1509 SurfPtr->n = (float) cos_alpha;
1510 if (fabs(SurfPtr->n) <= DBL_EPSILON) {
1511 SurfPtr->n = 0.0;
1512 }
1513
1514 #if 0
1515 if (Flags.warp == ON) {
1516 TileDisk1( warpfactor[i], SurfPtr );
1517 TileDisk2( warpfactor[i], SurfPtr );
1518 }
1519
1520 if (tilt_flag) {
1521 TileDisk1( TiltAngle, SurfPtr );
1522 TileDisk2( TiltAngle, SurfPtr );
1523 }
1524
1525 /* Calculate the layer vectors
1526 */
1527 CalculateLayerVecE1( SurfPtr );
1528 CalculateLayerVecE2( SurfPtr );
1529 #endif
1530
1531 SurfPtr->AreaType = TOP_SEGMENT;
1532
1533 /* Angle between x-Axis and Middle of the Element
1534 */
1535 SurfPtr->AreaAngle = j*360/n_phi;
1536 SurfPtr->MyNum = NumElem;
1537
1538 SurfPtr->rho = 1.00; /* ???*/
1539 SurfPtr->SpotNum = 0;
1540 SurfPtr->rad = r_j; /* Is needed for eclipse test */
1541
1542 SurfPtr->area = DSurfCor;
1543 Surf = Surf + DSurfCor;
1544
1545 /* --------------- INITIALIZE SURFACE TEMP ----------------- */
1546
1547 SurfPtr->temp = t_temp;
1548
1549 hs_phi = (phi < 0) ? (2.0 * PI + phi) : phi;
1550 hs_phi = fmod(hs_phi, 2.0 * PI);
1551
1552 hs_dist = DIST_CIRCLE(hs_phi, hs_long);
1553
1554 if ( (hs_dist <= hs_ext) && (r_j >= (Rout - hs_depth)) )
1555 {
1556 SurfPtr->temp = (hs_temp > t_temp) ? hs_temp : t_temp;
1557 SurfPtr->SpotNum = 0; /* FIXME 0 -> 1 */
1558 if (hs_puff > DBL_EPSILON)
1559 {
1560 double puffed_surf = puff_element(SurfPtr,
1561 (hs_dist/hs_ext),
1562 Rout,
1563 phi,
1564 h_i_old,
1565 h_i_new,
1566 DSurf /* to compute corrected area */);
1567 Surf -= DSurfCor;
1568 Surf += puffed_surf;
1569 }
1570 }
1571
1572 ++NumElem; /* increment surface element count */
1573 ++SurfPtr;
1574
1575 /* some elements of the structure are not filled at the moment */
1576
1577 } /* end of loop over phi */
1578
1579 r_i_old = r_i_new; /* end of one ring */
1580
1581 } /* end of loop over rings (n_r) */
1582
1583 /**********************************************************************/
1584 /* */
1585 /* ---------------------- Calculate outer side of disk ------------- */
1586 /* */
1587 /**********************************************************************/
1588
1589 /* area is outer side divided by rings/segments
1590 */
1591 DSurf = ((2.0*PI*Rout*Thick_out) / n_phi) / BinPtr->RingsOut;
1592
1593 /* vertical offset unit (half ring height)
1594 */
1595 offset_out = (Thick_out/BinPtr->RingsOut)/2.0;
1596
1597 /* ----------------- loop over rings (n_r) ---------------------- */
1598
1599 for (i = 0; i < BinPtr->RingsOut ; ++i) {
1600
1601 double zero_ring_height = ( Thick_out / 2.0 );
1602
1603 phi = 0.0; /* Correct ?????*/
1604
1605 /* Height of this ring. Start at top,
1606 * subtract half ring height, subtract (ring height * #ring)
1607 */
1608 ring_height = zero_ring_height
1609 - offset_out
1610 - i * (Thick_out/BinPtr->RingsOut);
1611
1612 if (Flags.debug[GEO]) {
1613 printf("OuterRing %03ld: %f\n", i, ring_height);
1614 }
1615
1616 /* check if enough memory */
1617 if ((NumElem+n_phi) >= (int) MaximumElements) {
1618
1619 if (Flags.debug[GEO])
1620 printf("Outer side of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements);
1621
1622 #ifdef _WITH_GTK
1623 if (Flags.interactive == OFF)
1624 nf_error(_(errmsg[10]));
1625 else {
1626 if (Flags.WantFit == OFF && Flags.WantMap == OFF)
1627 make_dialog(_(errmsg[10]));
1628 return(8);
1629 }
1630 #else
1631 nf_error(_(errmsg[10]));
1632 #endif
1633 }
1634
1635 /* dphi = i * ((2.0*PI/n_phi)/BinPtr->RingsOut); */
1636
1637 /* shift rings in phi by random amount to reduce numerical artifacts
1638 */
1639 #ifdef _WITH_OPENGL
1640 if (Flags.use_opengl == OFF)
1641 #endif
1642 dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0)));
1643
1644 if (t_flag < 2) /* isothermal */
1645 t_temp = T1;
1646 else
1647 t_temp = t_c * pow(Rout, -0.75);
1648
1649 /* ----------------- loop over phi -------------------------------- */
1650
1651 for (j = 0; j < n_phi; ++j) {
1652
1653 phi = ((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */
1654
1655 /* Shift rings in phi to avoid spikes in light curve.
1656 */
1657 phi += dphi;
1658
1659 sin_phi = sin(phi);
1660 cos_phi = cos(phi);
1661
1662 SurfPtr->RadiusIn = Rout; /* PR added 07.05.02 is needed for */
1663 SurfPtr->RadiusOut = Rout; /* eclipsing by disk verification */
1664
1665 SurfPtr->lambda = (float) Rout * cos_phi;
1666 SurfPtr->mu = (float) Rout * sin_phi;
1667 SurfPtr->nu = (float) ring_height;
1668
1669 SurfPtr->l = (float) cos_phi;
1670 if (fabs(SurfPtr->l) <= DBL_EPSILON){
1671 SurfPtr->l = 0.0;
1672 }
1673 SurfPtr->m = (float) sin_phi; // PR -sin -> sin 21.05.02
1674 if (fabs(SurfPtr->m) <= DBL_EPSILON){
1675 SurfPtr->m = 0.0;
1676 }
1677 SurfPtr->n = 0.0;
1678
1679 #if 0
1680 if (tilt_flag) {
1681 TileDisk1(TiltAngle,SurfPtr);
1682 TileDisk2(TiltAngle,SurfPtr);
1683 }
1684
1685 /* Calculate the layer vectors */
1686 CalculateLayerVecE1(SurfPtr);
1687 CalculateLayerVecE2(SurfPtr);
1688 #endif
1689
1690 SurfPtr->AreaType = OUT_RECTANGLE;
1691
1692 /* Angle between x-Axis and Middle of the Element
1693 */
1694 SurfPtr->AreaAngle = (j*360.0)/n_phi;
1695
1696 SurfPtr->AreaHeight = Thick_out/BinPtr->RingsOut;
1697 SurfPtr->AreaWidth = (2*PI*Rout)/n_phi;
1698 SurfPtr->MyNum = NumElem;
1699 SurfPtr->ring = i;
1700 SurfPtr->rad = Rout; /* Is needed for eclipse test */
1701
1702 SurfPtr->area = DSurf;
1703 Surf = Surf + DSurf;
1704
1705 /* --------------- INITIALIZE SURFACE TEMP ----------------- */
1706
1707 SurfPtr->temp = t_temp;
1708
1709 hs_phi = (phi < 0) ? (2.0 * PI + phi) : phi;
1710 hs_phi = fmod(hs_phi, 2.0 * PI);
1711
1712 hs_dist = DIST_CIRCLE(hs_phi, hs_long);
1713
1714 if ( hs_dist <= hs_ext )
1715 {
1716 SurfPtr->temp = (hs_temp > t_temp) ? hs_temp : t_temp;
1717 SurfPtr->SpotNum = 0; /* FIXME 0 -> 1 */
1718 if (hs_puff > DBL_EPSILON)
1719 {
1720 (void) puff_element(SurfPtr,
1721 (hs_dist/hs_ext),
1722 Rout,
1723 phi,
1724 zero_ring_height,
1725 ring_height,
1726 DSurf /* to compute corrected area */);
1727 }
1728 }
1729
1730 ++NumElem; /* increment surface element count */
1731 ++SurfPtr;
1732
1733 /* some elements of the structure are not filled at the moment */
1734
1735 } /* end of loop over phi */
1736 } /* end of loop over rings (n_r) */
1737
1738 /**********************************************************************/
1739 /* */
1740 /* ---------------------- Calculate bottom of disk ----------------- */
1741 /* */
1742 /**********************************************************************/
1743
1744 /* Here we loop from outer to inner border
1745 */
1746 r_outer_sqr = (r_i_old*r_i_old);
1747
1748 /* Area of an element. This can be lifted outside the loop, since
1749 * the progression of radii ensures equal areas for each element.
1750 */
1751 DSurf = (PI*P_0)/n_phi;
1752
1753 /* ----------------- loop over rings (n_r) ---------------------- */
1754
1755 for (i = 0; i < n_r; ++i) {
1756
1757 phi = 0.0;
1758
1759 /* check if enough memory */
1760 if ((NumElem+n_phi) >= (int) MaximumElements) {
1761
1762 if (Flags.debug[GEO])
1763 printf("Bottom of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements);
1764
1765 #ifdef _WITH_GTK
1766 if (Flags.interactive == OFF)
1767 nf_error(_(errmsg[10]));
1768 else {
1769 if (Flags.WantFit == OFF && Flags.WantMap == OFF)
1770 make_dialog(_(errmsg[10]));
1771 return(8);
1772 }
1773 #else
1774 nf_error(_(errmsg[10]));
1775 #endif
1776 }
1777
1778 r_i_new = sqrt(r_outer_sqr - (i+1)*P_0); /* radius of next circle */
1779 r_j = (r_i_new+r_i_old)/2.0; /* Center of the element */
1780
1781 /* warp_no = (int) n_r-i; */
1782
1783 ring_height = -1.0 * ((diskAdd + diskProp * pow(r_j, diskExp)) / 2.0);
1784
1785 /* angle of rise for this ring
1786 */
1787 h_i_old = -1.0 * ( (diskAdd + diskProp * pow(r_i_old, diskExp)) / 2.0);
1788 h_i_new = -1.0 * ( (diskAdd + diskProp * pow(r_i_new, diskExp)) / 2.0);
1789
1790 rise_angle = atan ( fabs(h_i_old - h_i_new) / fabs(r_i_old - r_i_new) );
1791
1792 /* Corrected surface area
1793 */
1794 DSurfCor = DSurf / cos ( rise_angle );
1795
1796 /* shift rings in phi by random amount to reduce numerical artifacts
1797 */
1798 #ifdef _WITH_OPENGL
1799 if (Flags.use_opengl == OFF)
1800 #endif
1801 dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0)));
1802
1803 /* Angle of rise from the middle to the edge of the Accretion Disk
1804 */
1805 if (t_flag == 0)
1806 alpha = (- 0.5 * diskProp);
1807 else
1808 alpha = -rise_angle /* (- 0.5 * diskProp) * diskExp * pow(r_j, (diskExp - 1.0)) */;
1809
1810 sin_alpha = sin(alpha);
1811 cos_alpha = cos(alpha);
1812
1813 if (t_flag < 2) /* isothermal */
1814 t_temp = T1;
1815 else
1816 t_temp = t_c * pow(r_j, -0.75);
1817
1818 /* ----------------- loop over phi -------------------------------- */
1819
1820 for (j = 0; j < n_phi; ++j) {
1821
1822 phi=((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */
1823
1824 phi += dphi;
1825
1826 sin_phi = sin(phi);
1827 cos_phi = cos(phi);
1828
1829 /* PR added 07.05.02 is needed for eclipsing
1830 * by disk verification
1831 * BEWARE!!! r_i_new und r_i_old muessen hier vertauscht sein
1832 * da die Schleife jetzt von Aussen nach Innen laeuft!!!!
1833 */
1834 SurfPtr->RadiusIn = r_i_new;
1835 SurfPtr->RadiusOut = r_i_old;
1836
1837 SurfPtr->lambda = (float) r_j * cos_phi;
1838 SurfPtr->mu = (float) r_j * sin_phi;
1839
1840 /* z negative (this is lower surface)
1841 */
1842 SurfPtr->nu = (float) ring_height;
1843
1844 SurfPtr->l = (float) (-cos_phi) * sin_alpha;
1845 if (fabs(SurfPtr->l) <= DBL_EPSILON) {
1846 SurfPtr->l = 0.0;
1847 }
1848 SurfPtr->m = (float) (-sin_phi) * sin_alpha;
1849 if (fabs(SurfPtr->m) <= DBL_EPSILON) {
1850 SurfPtr->m = 0.0;
1851 }
1852 SurfPtr->n = (float) -1.0 * (cos_alpha);
1853 if (fabs(SurfPtr->n) <= DBL_EPSILON) {
1854 SurfPtr->n = 0.0;
1855 }
1856
1857 #if 0
1858 if (Flags.warp == ON) {
1859 TileDisk1(warpfactor[warp_no],SurfPtr);
1860 TileDisk2(warpfactor[warp_no],SurfPtr);
1861 }
1862
1863 if (tilt_flag) {
1864 TileDisk1(TiltAngle,SurfPtr);
1865 TileDisk2(TiltAngle,SurfPtr);
1866 }
1867
1868 /* Calculate the layer vectors
1869 */
1870 CalculateLayerVecE1(SurfPtr);
1871 CalculateLayerVecE2(SurfPtr);
1872 #endif
1873
1874 SurfPtr->AreaType = BOTTOM_SEGMENT;
1875
1876 /* Angle between x-Axis and Middle of the Element
1877 */
1878 SurfPtr->AreaAngle = j*360/n_phi;
1879 SurfPtr->MyNum = NumElem;
1880 SurfPtr->ring = (n_r - i) - 1;
1881 SurfPtr->rad = r_j; /* Is needed for eclipse test */
1882 SurfPtr->area = DSurfCor;
1883 Surf = Surf + DSurfCor;
1884
1885 /* --------------- INITIALIZE SURFACE TEMP ----------------- */
1886
1887 SurfPtr->temp = t_temp;
1888
1889 hs_phi = (phi < 0) ? (2.0 * PI + phi) : phi;
1890 hs_phi = fmod(hs_phi, 2.0 * PI);
1891
1892 hs_dist = DIST_CIRCLE(hs_phi, hs_long);
1893
1894 if ( (hs_dist <= hs_ext) && (r_j >= (Rout - hs_depth) ))
1895 {
1896 SurfPtr->temp = (hs_temp > t_temp) ? hs_temp : t_temp;
1897 SurfPtr->SpotNum = 0; /* FIXME 0 -> 1 */
1898 if (hs_puff > DBL_EPSILON)
1899 {
1900 double puffed_surf = puff_element(SurfPtr,
1901 (hs_dist/hs_ext),
1902 Rout,
1903 phi,
1904 h_i_old,
1905 h_i_new,
1906 DSurf /* to compute corrected area */);
1907 Surf -= DSurfCor;
1908 Surf += puffed_surf;
1909 }
1910 }
1911
1912 ++NumElem; /* increment surface element count */
1913 ++SurfPtr;
1914
1915 /* some elements of the structure are not filled at the moment */
1916
1917 } /* end of loop over phi */
1918 r_i_old = r_i_new; /* end of one ring */
1919 } /* end of loop over rings (n_r) */
1920
1921
1922 /**********************************************************************/
1923 /* */
1924 /* ---------------------- Calculate inner side of disk ------------- */
1925 /* */
1926 /**********************************************************************/
1927
1928 /* area is inner side divided by rings/segments
1929 */
1930 DSurf = ((2.0*PI*Rin*Thick_in)/n_phi) / BinPtr->RingsIn;
1931
1932 /* vertical offset unit (half ring height)
1933 */
1934 offset_in = (Thick_in/BinPtr->RingsIn)/2.0;
1935
1936 /* ----------------- loop over rings (n_r) ---------------------- */
1937
1938 for (i = 0; i < BinPtr->RingsIn ; ++i) {
1939
1940 phi = 0.0; /* Correct ?????*/
1941
1942 /* Height of this ring. Start at bottom,
1943 * add half ring height, add (ring height * #ring)
1944 */
1945 ring_height = -1.0 * ( (diskAdd + diskProp * pow(Rin, diskExp)) / 2.0 )
1946 + offset_in
1947 + i * (Thick_in/BinPtr->RingsIn);
1948
1949 if (Flags.debug[GEO]) {
1950 printf("InnerRing %0ld: %f\n", i, ring_height);
1951 }
1952
1953 /* check if enough memory */
1954 if ((NumElem+n_phi) >= (int) MaximumElements) {
1955
1956 if (Flags.debug[GEO])
1957 printf("Inner side of disk: %ld >= %d\n", (NumElem+n_phi), (int) MaximumElements);
1958
1959 #ifdef _WITH_GTK
1960 if (Flags.interactive == OFF)
1961 nf_error(_(errmsg[10]));
1962 else {
1963 if (Flags.WantFit == OFF && Flags.WantMap == OFF)
1964 make_dialog(_(errmsg[10]));
1965 return(8);
1966 }
1967 #else
1968 nf_error(_(errmsg[10]));
1969 #endif
1970 }
1971
1972 /* dphi = i * ((2.0*PI/n_phi)/BinPtr->RingsIn); */
1973
1974 /* shift rings in phi by random amount to reduce numerical artifacts
1975 */
1976 #ifdef _WITH_OPENGL
1977 if (Flags.use_opengl == OFF)
1978 #endif
1979 dphi = (2.0*PI * (rand() / (RAND_MAX + 1.0)));
1980
1981 if (t_flag < 2) /* isothermal */
1982 t_temp = T1;
1983 else
1984 t_temp = t_c * pow(Rin, -0.75);
1985
1986
1987 /* ----------------- loop over phi -------------------------------- */
1988
1989 for (j = 0; j < n_phi; ++j) {
1990
1991 phi = ((2.0*PI/n_phi)*j)-Phase; /* MPR 23.05.2002 */
1992
1993 /* Shift rings in phi to avoid spikes in light curve.
1994 */
1995 phi += dphi;
1996
1997 sin_phi = sin(phi);
1998 cos_phi = cos(phi);
1999
2000 SurfPtr->RadiusIn = Rin; /* PR added 07.05.02 is needed for */
2001 SurfPtr->RadiusOut = Rin; /* eclipsing by disk verification */
2002
2003 SurfPtr->lambda = (float) Rin * cos_phi;
2004 SurfPtr->mu = (float) Rin * sin_phi;
2005 SurfPtr->nu = (float) ring_height;
2006
2007 /* Invert sign of cos(phi),/sin(phi), since the normal vector
2008 * points towards the z-axis here.
2009 */
2010 SurfPtr->l = (float) -1.0 * (cos_phi);
2011 if (fabs(SurfPtr->l) <= DBL_EPSILON) {
2012 SurfPtr->l = 0.0;
2013 }
2014 SurfPtr->m = (float) -1.0 * (sin_phi);
2015 if (fabs(SurfPtr->m) <= DBL_EPSILON) {
2016 SurfPtr->m = 0.0;
2017 }
2018 SurfPtr->n = 0.0;
2019
2020 #if 0
2021 if (tilt_flag) {
2022 TileDisk1(TiltAngle,SurfPtr);
2023 TileDisk2(TiltAngle,SurfPtr);
2024 }
2025
2026 /* Calculate the layer vectors */
2027 CalculateLayerVecE1(SurfPtr);
2028 CalculateLayerVecE2(SurfPtr);
2029 #endif
2030
2031 SurfPtr->AreaType = IN_RECTANGLE;
2032
2033 /* Angle between x-Axis and Middle of the Element
2034 */
2035 SurfPtr->AreaAngle = j*360/n_phi;
2036 SurfPtr->AreaHeight = Thick_in/BinPtr->RingsIn;
2037 SurfPtr->AreaWidth = (2*PI*Rin)/n_phi;
2038 SurfPtr->MyNum = NumElem;
2039 SurfPtr->ring = i;
2040 SurfPtr->rad = Rin; /* Is needed for eclipse test */
2041
2042 SurfPtr->area = DSurf;
2043 Surf = Surf + DSurf;
2044
2045 /* --------------- INITIALIZE SURFACE TEMP ----------------- */
2046
2047 SurfPtr->temp = t_temp;
2048
2049 hs_phi = (phi < 0) ? (2.0 * PI + phi) : phi;
2050 hs_phi = fmod(hs_phi, 2.0 * PI);
2051
2052 hs_dist = DIST_CIRCLE(hs_phi, hs_long);
2053
2054 if ( (hs_dist <= hs_ext) &&
2055 (Rin >= (Rout - hs_depth) ))
2056 {
2057 SurfPtr->temp = (hs_temp > t_temp) ? hs_temp : t_temp;
2058 SurfPtr->SpotNum = 0; /* FIXME 0 -> 1 */
2059 }
2060
2061 ++NumElem; /* increment surface element count */
2062 ++SurfPtr;
2063
2064 /* some elements of the structure are not filled at the moment */
2065
2066 } /* end of loop over phi */
2067 } /* end of loop over rings (n_r) */
2068
2069 BinPtr->NumElem = NumElem; /* Number of elements*/
2070
2071 if (Flags.debug[VERBOSE] == ON) {
2072 printf("\n");
2073 printf(_(" Divide Disk %3d : %6ld Surface elements\n"),
2074 Comp+1, BinPtr->NumElem);
2075 }
2076
2077 BinPtr->Surface = Surf; /*Added 29-04-2002, EG/MPR */
2078
2079 /* ---------------------------- mean temperature ------------------- */
2080
2081 Weight = 0.0; TempMean = 0.0;
2082 SurfPtr = Surface[Comp];
2083 for (i = 0; i < BinPtr->NumElem; ++i) {
2084 Area = SurfPtr->area;
2085 TempMean = TempMean + SurfPtr->temp * Area;
2086 Weight = Weight + Area;
2087 /*------- Sollte hier nicht irgendwo ein SurfPtr++ stehen????? */
2088 /* eigentlich schon ... RW Tue Sep 11 12:54:20 CEST 2001 */
2089 ++SurfPtr;
2090 }
2091
2092 BinPtr->TempMean = TempMean/Weight;
2093
2094 /**********************************************************************/
2095 /* */
2096 /* Flux intensity computation */
2097 /* */
2098 /**********************************************************************/
2099
2100 if (Flags.blackbody == ON)
2101 {
2102 BlackbodyFlux(Disk);
2103 }
2104 else
2105 {
2106 Binary[Disk].log_g = 3.5;
2107 ModelFlux(Disk);
2108 }
2109
2110 /* ---------------- output --------------------- */
2111
2112 if (Flags.debug[VERBOSE] == ON) {
2113 /* printf("\n Star: %5d\n", (Comp+1)); */
2114 printf(_(" Surface Area, Mean Gravity (dimensionless): %8.6f %8.6f\n"),
2115 BinPtr->Surface, BinPtr->Gravity);
2116 printf(_(" Temperature is %8.6f\n"), BinPtr->Temperature);
2117
2118 printf("\n");
2119 }
2120
2121 return(0);
2122 }
2123 #endif
2124
2125
2126
2127
2128