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 #include "config.h" 21 22 /* This is actually nowhere used, but causes a 'redefined' */ 23 /* error due to a bug in the JPEG lib headers. */ 24 #undef HAVE_STDLIB_H 25 26 /* Fourth LD law is currently for experimental purposes only */ 27 #define LD_LAW_MAX 3 28 29 /* --------------------------------------------------------- */ 30 /* */ 31 /* Default values that may be changed */ 32 /* */ 33 /* --------------------------------------------------------- */ 34 35 /* default number of steps per pi (in Eta) */ 36 /* pi = 3.14 = 180 degree */ 37 /* steps in Phi are calculated as 10+2*STEPS_PER_PI/sin(Eta) */ 38 /* -- and corresponding maximum number of surface elements */ 39 /* */ 40 41 #ifdef HIGH_PRECISION 42 # define STEPS_PER_PI 196 43 # define MAXELEMENTS 65536 44 #else 45 # define STEPS_PER_PI 64 46 # ifdef HAVE_DISK 47 # define MAXELEMENTS 65536 48 # else 49 # define MAXELEMENTS 8192 50 # endif 51 #endif 52 53 54 /* window size if Gnuplot */ 55 /* for pgplot, width will be used to scale display */ 56 /* unless environment variable PGPLOT_XW_WIDTH is set */ 57 /* */ 58 #define GNU_GEOMETRY "550x424+300+20" 59 60 /* maximum number of steps for a full orbit */ 61 /* */ 62 #define PHASESTEPS 8192 63 64 /* minimum number of steps for a full orbit */ 65 /* */ 66 #ifdef HIGH_PRECISION 67 # define DEF_STEPS 360 68 #else 69 # define DEF_STEPS 80 70 #endif 71 72 #define FILTERNAME_SIZE 16 73 74 #define LOW_STEPS 8 75 76 /* maximum number of observations in each passband/filter */ 77 /* */ 78 #define MAXOBS 8192 79 80 /* maximum number of spots per component */ 81 /* */ 82 #define N_SPOT 10 83 84 /* output file for lightcurve */ 85 /* */ 86 #define OUT_FILE "NightfallCurve.dat" 87 88 /* output file for fit results */ 89 /* */ 90 #define OUT_SIMPLEX "NightfallFitted.dat" 91 92 /* output file for chi square map */ 93 /* */ 94 #define OUT_CHIMAP "NightfallChiMap.dat" 95 96 /* Grid size (square) for Chi Square Map */ 97 /* */ 98 #define _CHI_SCANS_ 16 99 100 /* default wavelength for line profile (V band) */ 101 /* */ 102 #define LAMDAZERO_DEFAULT 550.0 103 104 #ifdef HIGH_PRECISION 105 /* array for line profile (multiple of two) */ 106 #define PROFILE_ARRAY 1024 107 /* resolution of profile (delta_lambda/lambda) */ 108 #define PROFILE_RES 160000.0 109 #else 110 /* array for line profile (multiple of two) */ 111 #define PROFILE_ARRAY 128 112 /* resolution of profile (delta_lambda/lambda) */ 113 #define PROFILE_RES 20000.0 114 #endif 115 116 /* switch model atmospheres PHOENIX -> ATLAS */ 117 #define SWITCH_TEMP 9000.0 118 119 /* --------------------------------------------------------- */ 120 /* */ 121 /* nothing to change below this line */ 122 /* (unless you know what you do ... ) */ 123 /* */ 124 /* --------------------------------------------------------- */ 125 126 /* Maximum Number of Iterations for RootFind/MinFind */ 127 #define MAXITER 100 128 129 /* Maximum linesize for config file */ 130 #define MAX_CFG_INLINE 255 131 132 /* Floating-Point precision */ 133 #define EPS 3.0e-8 134 135 /* number of passbands */ 136 #define NUM_MAG 12 137 138 /* maximum iteration for downhill simplex */ 139 #define MAX_ITER_SIMPLEX 200 140 141 /* default tolerance for downhill simplex */ 142 #define SIMPLEX_FIT_TOLERANCE 0.01 143 144 /* the maximum dimension of the parameter space */ 145 /* -- Careful: some calls to array indices are */ 146 /* hardcoded (LightSimplex.c), decreasing SDIM */ 147 /* can lead to disaster */ 148 149 #define SDIM 70 150 151 /* Size of Flags.simplex */ 152 #define MAX_SPX 128 153 154 /************************************************************************/ 155 /* */ 156 /* TYPEDEFS */ 157 /* */ 158 /************************************************************************/ 159 160 #define WARN 0 161 #define BUSY 1 162 #define VERBOSE 2 163 #define GEO 3 164 #define OPTIMIZE 4 165 #define MINFIND 5 166 #define CHISQR 6 167 #define INFILE 7 168 #define SURFACE 8 169 #define TERSE 9 170 #define PARALLEL 9 171 /* Array bound for above flags */ 172 #define MAX_FVAL 10 173 174 #define NF_NATIVE_FILE 0 175 176 /* ----------------- A 2D Point ------------------------------- */ 177 178 typedef struct Point2Struct { /* 2D Point */ 179 double x, y; 180 } Point2; 181 typedef Point2 Vector2; 182 183 /* ----------------- Output Table ------------------------------- */ 184 185 typedef struct PhotoOutStruct { 186 float Phase; /* phase of output datum */ 187 double Mag[NUM_MAG]; /* magnitudes of output data */ 188 float RV[3]; /* radial velocities */ 189 /* 2->3 components 2002/04/09, MPR */ 190 float D_RV[3]; /* disturbance of Kepler RV */ 191 /* 2->3 components 2002/04/09, MPR */ 192 } PhotoOut; 193 194 /* ----------------- Data Point ------------------------------- */ 195 196 typedef struct PhotoStruct { 197 double time; /* time of observation */ 198 float mag; /* value of observation */ 199 float weight; /* weight of observation */ 200 float phi; /* phase of observation */ 201 float residual; /* residual of observation */ 202 } Photo3; 203 204 /* ----------------- Spot Parameters ------------------------------- */ 205 206 typedef struct SpotStruct { 207 double nowlongitude; /* actual value */ 208 double longitude; /* latitude of spot */ 209 double latitude; /* longitude of spot */ 210 double radius; /* angular radius */ 211 double dimfactor; /* dimming (= T.new/T.old) */ 212 } SpotLight; 213 214 /* ----------------- Surface Element ------------------------------- */ 215 216 typedef struct SurfaceStruct { 217 float rad; /* radius */ 218 float eta; /* angle to x-axis */ 219 float phi; /* angle in y-z-plane */ 220 float lambda; /* holds rad*cos(eta) = x */ 221 float mu; /* holds rad*sin(eta)cos(phi) = y */ 222 float nu; /* holds rad*sin(eta)sin(phi) = z */ 223 float l; /* surface normal x */ 224 float m; /* surface normal y */ 225 float n; /* surface normal z */ 226 float rho; /* distance to other stars centre */ 227 #ifdef HIGH_PRECISION 228 double area; /* area of element */ 229 #else 230 float area; /* area of element */ 231 #endif 232 float grav; /* dimensionless gravity at element*/ 233 double temp; /* temperature of element */ 234 double orig_temp; /* temperature of element */ 235 double f_[NUM_MAG]; /* intensity in various passbands */ 236 long MyNum; /* number */ 237 int ring; /* ring this element lives in */ 238 struct SurfaceStruct *phi_ptr; /* pointer to last in phi */ 239 struct SurfaceStruct *phi1_ptr; /* pointer to next in phi */ 240 struct SurfaceStruct *next_ptr; /* pointer to next in eta */ 241 struct SurfaceStruct *last_ptr; /* pointer to last in eta */ 242 float LOS_Pot; /* Potential min. on LOS */ 243 float visibility; /* (fractional) visibility */ 244 float CosGamma; /* angle LOS-surface normal */ 245 int EclipseFlag; /* eclipsed = ON/OFF */ 246 int SpotNum; /* number of overlapping spots */ 247 double SumDimFactor; /* average dimming factor */ 248 float Velocity; /* projected velocity */ 249 250 /* Parameter needed for eclipse testing */ 251 int AreaType; /* AreaType: 252 (TOP_SEGMENT/BOTTOM_SEGMENT/ 253 OUT_RECTANGLE/IN_RECTANGLE) */ 254 double RadiusIn; /* Inner radius of area element PR */ 255 double RadiusOut; /* Outer radius of area element PR */ 256 double AreaAngle; /* Angle of area element between 257 x-axis and middle of the element*/ 258 double AreaHeight; 259 double AreaWidth; 260 double VecE1[3]; /* Element vector 1 of test layer */ 261 double VecE2[3]; /* Element vector 2 of test layer */ 262 } SurfaceElement; 263 264 /* ----------------- Disk Component -------------------------------- */ 265 /* disk parameters MK 14.03.2001 */ 266 267 /* typedef struct DiskStruct { */ 268 /* double Incl; */ /* disk inclination */ 269 /* */ /* in degrees */ 270 /* double Warp; */ /* warp angle of the disk */ 271 /* */ /* in degrees */ 272 /* double Thick; */ /* thickness of disk */ 273 /* */ /* in units of R */ 274 /* double HR; */ /* H over R for the disk */ 275 /* double Rout; */ /* outer radius of the disk */ 276 /* */ /* in units of seperation */ 277 /* double Rin; */ /* inner radius of the disk */ 278 /* */ /* in units of seperation */ 279 /* } DiskComponent; */ 280 281 /* ----------------- Binary Component ----------------------------- */ 282 283 typedef struct ComponentStruct { 284 double Mq; /* mass ratio */ 285 double RocheFill; /* Roche Lobe filling Factor */ 286 double RocheStore; /* Roche Lobe filling Factor */ 287 double RLag1; /* Lagrange One */ 288 double RLag2; /* Lagrange Two */ 289 double RLag3; /* Lagrange Three */ 290 double RCLag1; /* Lagrange One Potential */ 291 double RCLag2; /* Lagrange Two Potential */ 292 double RCLag3; /* Lagrange Three Potential */ 293 double CosPhiLagrange; /* Lim. angle for eclipse */ 294 double RZatL1; /* Throat Rad perpend. L1 */ 295 double RYatL1; /* Throat Rad perpend. L1 */ 296 double LimRad; /* Radius at Throat */ 297 double LimEta; /* Eta at Throat */ 298 double Fratio; /* Nonsync Rotation */ 299 double RXCrit; /* Max. Radius on X axis */ 300 double RochePotCrit; /* Potential for Max Lobe */ 301 double RochePot; /* Roche Potential-Star */ 302 double RCrit; /* Radius for Max Lobe */ 303 double Radius; /* Polar Radius Star */ 304 double Surface; /* Surface Area */ 305 double Gravity; /* MeanGravity */ 306 double Temperature; /* Effective Temperature */ 307 double Albedo; /* Bolometric Albedo */ 308 double GravDarkCoeff; /* Gravitational Darkening */ 309 double DarkeningGrav; /* Denominator for dark. law */ 310 long NumElem; /* Number of Surface Elements */ 311 long NumPlus; /* plus the Throat ... */ 312 long N_PhiStep[STEPS_PER_PI+10]; 313 /* steps in phi at eta */ 314 int N_Rings; /* number of rings (in eta) */ 315 /* steps in phi at eta */ 316 double Xp; /* Smallest enclosing radius */ 317 double Volume; /* volume of star */ 318 double ScaledVolume; /* scaled volume */ 319 double RadMean; /* mean radius */ 320 double TempMean; /* mean temperature */ 321 double K; /* half-amplitude of radial v */ 322 double log_g; /* log_g for model atmosphere */ 323 324 /* -------------------- disk parameters -------------------------- */ 325 326 double Tilt; /* disk inclination */ 327 /* in degrees */ 328 double Warp; /* warp angle of the disk */ 329 /* in degrees */ 330 double Thick; /* thickness of disk */ 331 /* in units of R */ 332 double HR; /* H over R for the disk */ 333 double Rout; /* outer radius of the disk */ 334 /* in units of seperation */ 335 double Rin; /* inner radius of the disk */ 336 /* in units of seperation */ 337 int RingsOut; /* Number of rings on the 338 outer side of the disk */ 339 int RingsIn; /* Number of Rings on the * 340 inner side of the disk */ 341 int Segments; /* Number of Segments in Phi 342 direction */ 343 344 /* -------------------- hot spot parameters ---------------------- */ 345 346 double tempHS; /* Hot Spot temperature */ 347 double longitudeHS; /* Hot Spot longitude */ 348 double longitude_intern;/* internal value */ 349 double extentHS; /* Hot Spot extent */ 350 double depthHS; /* Hot Spot depth */ 351 double puffHS; /* Hot Spot puff */ 352 353 /* -------------------- warp parameters -------------------------- */ 354 355 double maxWarp; /* warp maximum */ 356 double longitudeWarp; /* warp longitude */ 357 double widthWarp; /* warp width */ 358 double centreWarp; /* warp centre */ 359 360 } BinaryComponent; 361 362 /* ----------------- Flags --------------------------------------- */ 363 /* added disk related flags MK 13.05.01 */ 364 365 typedef struct FlagsStruct { 366 unsigned char debug[MAX_FVAL]; /* debug flags */ 367 int disk; /* want disk geometry (ON/OFF) */ 368 int warp; /* warped disk (ON/OFF) */ 369 int plot; /* want lightcurve (ON/OFF) */ 370 int visualize; /* want visualize (ON/OFF) */ 371 int animate; /* want animated (ON/OFF) */ 372 #ifdef _WITH_OPENGL 373 int axes; /* want x-,y-, z-axis (ON/OFF) */ 374 int movie; /* want jpeg movie (ON/OFF) */ 375 int texture; /* want textures (ON/OFF) */ 376 int textype; /* what kind of textures (IMAGE,TEMP)*/ 377 int frame; /* number of jpeg frame */ 378 int wireframe; /* want wire frame or solid rendering*/ 379 int points; /* use POINTS instead of QUAD_STRIP */ 380 int labels; /* want labels (ON/OFF) */ 381 float GLdistancez; /* distance of observer in */ 382 /* z-direction */ 383 float GLdistancey; /* distance of observer in */ 384 /* y-direction */ 385 float GLdistancex; /* distance of observer in */ 386 /* x-direction */ 387 int use_opengl; /* really use it ? (ON/OFF) */ 388 #endif 389 int interactive; /* want interactive (ON/OFF) */ 390 int edu; /* educational */ 391 int fill; /* overfill roche lobe (ON/OFF) */ 392 int lineprofile; /* want line profile (ON/OFF) */ 393 int smap; /* want surface map (ON/OFF) */ 394 int monochrome; /* want monochrome wavelen (ON/OFF) */ 395 int limb; /* limb darkening method */ 396 int reflect; /* reflection treatment */ 397 int fractional; /* fractional visibility (ON/OFF) */ 398 int blackbody; /* blackbody approximation (ON/OFF) */ 399 int InEclipse1; /* we are in eclipse w/ Primary */ 400 int InEclipse2; /* we are in eclipse w/ Secondary */ 401 int asynchron1; /* Primary rotates asynchron */ 402 int asynchron2; /* Secondary rotates asynchron */ 403 int Spots1; /* num. of spots on primary */ 404 int Spots2; /* num. of spots on secondary */ 405 int elliptic; /* is eccentric orbit */ 406 int first_pass; /* first pass in main_loop() */ 407 long Passbands[NUM_MAG+2]; /* number of data read in passband N */ 408 int PlotBand; /* which passband to plot ? */ 409 unsigned char simplex[MAX_SPX]; /* parameter set to fit */ 410 int WantFit; /* want to fit something (ON/OFF) */ 411 float SimplexTol; /* how good is the fit ? */ 412 float Step[2]; /* Stepsizes for Map */ 413 int WantMap; /* we want to map chi square */ 414 int anneal; /* do simulated annealing */ 415 int eps; /* EPS plot */ 416 int IsComputed; /* Lightcurve is computed */ 417 int ProComputed; /* Profile is computed */ 418 int parseCL; /* do we parse the command line */ 419 int plotOpened; /* Gnuplot window opened */ 420 int GN_exitconfirm; /* whether confirm exit */ 421 int GN_rempos; /* whether remember position */ 422 int GN_tooltips; /* whether have tooltips */ 423 int GL; /* does the display support OpenGL */ 424 char ConfFile[256]; /* configuration file */ 425 int tdisk; /* disk temperature distribution */ 426 } FlagsHandle; 427 428 /* ----------------- The Orbit ------------------------------------- */ 429 430 typedef struct OrbitStruct { 431 char Name[MAX_CFG_INLINE+1]; /* system name */ 432 double Inclination; /* orbit inclination */ 433 double Omega; /* length of periastron */ 434 double OmegaZero; /* passage of periastron */ 435 double MQuad; /* start (quadrature) */ 436 double Excentricity; /* excentricity of orbit */ 437 double Nu; /* true anomaly */ 438 double M; /* mean anomaly */ 439 double E; /* excentric anomaly */ 440 double MinDist; /* distance at periastron */ 441 double Phase; /* the phase angle */ 442 double Dist; /* the distance */ 443 double TruePeriod; /* Period in seconds */ 444 double TrueMass; /* Total mass in kg */ 445 double TrueDistance; /* Periastron Dist in meter */ 446 double LambdaZero; /* Wavelength for Profile */ 447 int PhaseIndex; /* Where do we stand ? */ 448 double Third[NUM_MAG]; /* Third Light */ 449 /* Contact[0] gets zeroed somewhere... valgrind can't find it :( */ 450 double Contact[9]; /* Contact phases */ 451 } OrbitType; 452 453 #define PERIDIST (Orbit.MinDist*Orbit.TrueDistance) 454 455 /* ----------------- Data Files ------------------------------------ */ 456 457 #define MAX_INFILE_PATH 4096 458 459 typedef struct FileStruct { 460 char DataFile[MAX_INFILE_PATH]; 461 int DataFormat; 462 /*@null@*/ struct FileStruct *nextFile; 463 } FileType; 464 465 /*********************************/ 466 /* one-argument macros */ 467 /*********************************/ 468 469 470 /* round a to nearest integer */ 471 #define ROUND(a) (((a) >= FLT_EPSILON) ? (int)((a)+0.5) : -(int)(0.5-(a))) 472 473 /* square a */ 474 extern double GArgSquare; 475 # define SQR(a) \ 476 ( fabs(GArgSquare = (a)) <= FLT_EPSILON ? 0.0 : (GArgSquare*GArgSquare) ) 477 478 # define SQR_DIRECT(a) \ 479 ( (a) <= FLT_EPSILON ? 0.0 : ((a)*(a)) ) 480 481 482 /*********************************/ 483 /* error handling */ 484 /*********************************/ 485 486 /* shout if something goes wrong */ 487 #define WARNING(msg) \ 488 {if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON) \ 489 fprintf(stderr, "**Warning**: %s \n", msg);} 490 491 #define FATAL(msg) \ 492 {fprintf(stderr, "**ERROR**: %s \n", msg); exit(EXIT_FAILURE);} 493 494 /* input error */ 495 #define INERR(msg, input) \ 496 {fprintf(stderr,"** Input Error **: %s \n (offending item: %s)\n .. exit to system\n", \ 497 msg,input); \ 498 exit(EXIT_FAILURE);} 499 500 /* input error */ 501 #define ERRLIM(msg, low, high, input) \ 502 {fprintf(stderr,"** Input Error **: %s \n (limits: %f..%f, offending item: %s)\n .. exit to system\n", \ 503 msg,low, high, input); \ 504 exit(EXIT_FAILURE);} 505 506 /* input warning */ 507 #define WARNERR(msg, input) \ 508 {fprintf(stderr, "** Input Warning **: %s \n (offending item: %s)\n", msg, input);} 509 510 /* input warning */ 511 #define WARNLIM(msg, low, high, input) \ 512 {fprintf(stderr, "** Input Warning **: %s \n (limits: %f..%f, offending item: %s)\n", msg, low, high, input);} 513 514 515 /*********************************/ 516 /* two-argument macros */ 517 /*********************************/ 518 519 /* find minimum of a and b */ 520 #undef MIN 521 #define MIN(a,b) ( ( (a)-(b) <= FLT_EPSILON ) ? (a) : (b)) 522 523 /* find minimum of a and b */ 524 #undef IMIN 525 #define IMIN(a,b) ( ( (a)<(b) ) ? (a) : (b)) 526 527 /* find maximum of a and b */ 528 #undef MAX 529 #define MAX(a,b) ( ( (a)-(b) >= FLT_EPSILON ) ? (a) : (b)) 530 531 /* find maximum of a and b */ 532 #undef IMAX 533 #define IMAX(a,b) ( ( (a)>(b) ) ? (a) : (b)) 534 535 /* Ensures that x is between the limits set by low and high */ 536 #undef CLAMP 537 #define CLAMP(x, low, high) \ 538 (((x) > (high)) ? (high) : (((x) < (low)) ? (low) : (x))) 539 540 541 /*********************************/ 542 /* useful constants */ 543 /*********************************/ 544 545 #ifndef PI 546 #define PI 3.14159265358979323846 547 #endif 548 #ifndef LN10 549 #define LN10 2.30258509299404568402 /* ln 10; 10eX = exp(ln10*X) */ 550 #endif 551 #ifndef DTOR 552 #define DTOR 0.01745329251994329547 /* radian = degree * DTOR */ 553 #endif 554 #ifndef RTOD 555 #define RTOD 57.2957795130823228646 /* degree = radian * RTOD */ 556 #endif 557 558 #ifndef EULER 559 #define EULER 2.71828182845904523536 560 #endif 561 562 #define MCMC 2 563 #define ON 1 564 #define OFF 0 565 566 /* Sun radius in meter 567 */ 568 #define SOL_RADIUS 6.960E8 569 570 /* Surface identifier for the disk */ 571 572 #define TOP_SEGMENT 1 573 #define OUT_RECTANGLE 2 574 #define BOTTOM_SEGMENT 3 575 #define IN_RECTANGLE 4 576 577 #define EclipsedByPrimary 1 578 #define EclipsedBySecondary 2 579 #define EclipsedByDisk 3 580 581 /* illuminated/non-illuminated part of stellar surface */ 582 583 #define BACKSIDE 1 584 #define FRONTSIDE (-1) 585 586 587 /*****************************************************/ 588 /* */ 589 /* PROTOTYPES */ 590 /* */ 591 /*****************************************************/ 592 593 594 /* ---------------- the gui --------------- */ 595 596 #ifdef _WITH_GTK 597 /* main GTK program */ 598 int the_gui(int argc, char *argv[]); 599 #endif 600 601 /* ------- temperature distribution --------------- */ 602 603 /* compute detailed reflection */ 604 void LightReflect(); 605 606 /* simple reflection treatment */ 607 void SimpleReflect(int Comp); 608 609 /* Compute Gravitation Darkening, Albedo */ 610 void ComputeGravDark(); 611 612 /* put spots on surface */ 613 void MakeSpots(int Component, int phasestep); 614 615 /* ------- flux computation ----------------------- */ 616 617 /* Get a set of monochromatic wavelenghts from the */ 618 /* environment */ 619 int DefineMonoWavelengths (); 620 621 /* compute temperature-dependent effective */ 622 /* wavelengths for blackbody */ 623 void EffectiveWavelength(int Comp); 624 625 /* Compute limb darkening coeff - filter dependedent */ 626 void LightLimbDark(int Comp); 627 628 /* flux in blackbody approximation */ 629 void BlackbodyFlux(int Comp); 630 631 /* flux with atmospheric model */ 632 void ModelFlux(int Comp); 633 634 /* limits for Kurucz models */ 635 float ModelLimits(float log_g); 636 637 /* flux computing - filter dependent */ 638 void LightFlux(int Comp, int Phase_index); 639 640 /* normalize the output flux */ 641 void NormalizeFlux(); 642 643 /* add the doppler boost */ 644 void DopplerBoostFlux(); 645 646 /* Initialize Output flux tables - filter dependent */ 647 void InitOutFlux(); 648 649 /* normalize the output flux */ 650 void SetNormalPhase (int band, float phase); 651 652 /* reset normalization point */ 653 void ResetNormalPhaseOne (int band); 654 655 /* reset all normalization points */ 656 void ResetNormalPhase (); 657 658 /* copy flux to L90 array */ 659 void LightLuminosity(int component, int phaseindex, int n_phase); 660 661 /* ------ simplex fit algorithm ------------------- */ 662 663 /* the downhill simplex code */ 664 int Simplex(); 665 666 /* morph Simplex */ 667 double SimplexFlow(double Factor, double M[][SDIM], 668 double Y[], double M_Sum[], int Worst, 669 int FitDim, int *ErrCode); 670 671 /* Initialize the ranges */ 672 void SimplexInitRanges(); 673 674 /* Initialize the Simplex */ 675 int SimplexInit(double Y[], double M[][SDIM]); 676 677 /* set the global variables */ 678 void SimplexSetVariables(double x[]); 679 680 /* check the Vertex */ 681 int SimplexCheckVertex(double X[]); 682 683 /* print fit result the global variables */ 684 void SimplexPrintVariables(double Chi, 685 double VarChiSquare); 686 687 /* simulated annealing */ 688 float SaAnneal(double tscale, double t0, 689 double X[], double *Y_Best, 690 double BestEver[], double *BestChi); 691 692 /* Markov Chain Monte Carlo */ 693 int MCMC_go(double X[], double *Y_Best, 694 double BestEver[], double *BestChi); 695 696 /* Map Chi Square */ 697 int ChiMap(); 698 699 /* ------ main ----------------------------------- */ 700 701 /* initialize flags and variables */ 702 void InitFlags(); 703 704 /* Print Usage Info */ 705 void UsageInfo(); 706 707 /* Get The Command Line Arguments */ 708 void GetArguments(int *argc, char *argv[]); 709 710 /* Here we set up the Stars, do the Lightcurve, */ 711 /* and return the Merit */ 712 int MainLoop(/*@out@*/ double *Merit); 713 714 /* radial velocity curve */ 715 void RadVel(int Phase_Index); 716 717 /* ------ distributed ---------------------------- */ 718 719 /* Return 1 if host file exist, else 0 */ 720 int have_parallel_file (); 721 722 /* Reset everything */ 723 void DistReset(); 724 725 /* Execute */ 726 int DistExecLoop(double * merit); 727 728 /* Collect results */ 729 int DistCollectMerits(double * merit); 730 731 /* ------ statistics ----------------------------- */ 732 733 /* compute the merit function - we use ChiSquare */ 734 /* globally over all input data */ 735 double MeritFunction (); 736 737 /* perform runs test on fit residuals, return number */ 738 /* of runs, upper and lower 5 percent limit */ 739 void Runs(Photo3 *data, long NumData, int *Runs, 740 int *UpLim, int *LowLim); 741 742 /* -------- Input/Output --------------------------- */ 743 744 /* wrapper function to do printing */ 745 void my_cpgend(); 746 747 /* Get data file path */ 748 char * data_data_fls(); 749 750 /* Get cfg file path */ 751 char * data_cfg_fls(); 752 753 /* Get doc file path */ 754 char * data_doc_fls(); 755 756 /* Get pixmap file path */ 757 char * data_pix_fls(); 758 759 /* parse the config file */ 760 void ParseConfig(const char *InputFile, int *numarg); 761 762 /* write header of output file */ 763 void OutfileHeader(FILE *file_out); 764 765 /* read a single line from an open file */ 766 int LireLigne(char *s, int lim, FILE *fpe); 767 768 /* read from input data file */ 769 void Read(const char InputFile[]); 770 771 /* visualize geometry */ 772 void PlotGeometry(); 773 774 /* Write Output flux tables - filter dependent */ 775 void WriteOutput(); 776 777 /* Plot Output flux tables - filter dependent */ 778 void PlotOutput(); 779 780 /* animated plot */ 781 void Animate(int j); 782 783 /* set PGPLOT window size */ 784 void SetPgplotWindow(); 785 786 /* --------- Eclipse verification, Visibility ----- */ 787 788 /* attach the throat of star a to star b */ 789 void LightCopyThroat(); 790 791 /* copy the throat back */ 792 void LightCopyThroatBack(); 793 794 /* Calculate various quantities for simple */ 795 /* geometric tests */ 796 int LightSetupTests(); 797 798 /* Eclipse Verification for stars */ 799 void LightCurve(int Phase_index); 800 801 /* Eclipse Verification for disk */ 802 void LightCurveDisk(int j); 803 void eclipsedbydisk(double lzero2, double mzero2, double nzero2, long Diskelements, int Comp, SurfaceElement *SurfPtrD, int *eclipsed); 804 805 /* Fractional Visibility */ 806 void LightFractionalPot(int Component); 807 808 /* -------- Functions ----------------------------- */ 809 810 /* find minimum of potential along LOS */ 811 int MinFinder(double *q, double *f, 812 double *t1, double *t3, 813 double *l0, double *m0, double *n0, 814 double *x0, double *y0, double *z0, 815 double *tmin, /*@out@*/ double *Cmin); 816 817 /* find root of Function in [low,high] w/ tolerance */ 818 double RootFind(double (*Function)(double), 819 const double Low, const double Up, 820 const double Tolerance, const char *caller, int *Err); 821 822 double RootBisect(double (*Function)(double), 823 const double Low, const double Up, 824 const double Tolerance, const char *caller, int *Err); 825 826 827 /* Zero at LagrangeOne */ 828 double LagrangeOne(double x); 829 830 /* Zero at LagrangeTwo */ 831 double LagrangeTwo(double x); 832 833 /* Zero at LagrangeThree */ 834 double LagrangeThree(double x); 835 836 /* Potential in X-Y Plane */ 837 double RocheXYPrimary(double x, double y, double q, 838 double f); 839 840 /* Potential */ 841 double RocheSurface(double x); 842 843 /* Potential */ 844 double RocheYPerpendL1(double y); 845 846 /* Potential */ 847 double RocheZPerpendL1(double z); 848 849 /* Potential */ 850 double RochePerpendSecond(double z); 851 852 /* find the phase of periastron */ 853 void FindPeriastron(); 854 855 /* solve kepler equation */ 856 void SolveKepler(int j); 857 858 /* shift phase to correct interval */ 859 void PhaseShift(); 860 861 /* analytic polar radius */ 862 double PolRad(double q, double r); 863 864 /* returns (volume-scaledvolume) as function of r */ 865 /* zero if r = RadMean(ScaledVolume) */ 866 double VolOfR(double r); 867 868 /* copies PhotoPoint b to a */ 869 void V2CopyPhoto(/*@out@*/ Photo3 *a, Photo3 *b); 870 871 /* Sort in ascending order in phi w/ */ 872 /* Heapsort Algorithm */ 873 /* D.E. Knuth, Sorting and Searching, 1973 */ 874 /* contrary to Quicksort, this is guaranteed to be */ 875 /* always of order NlogN, even in worst case */ 876 void SortPhotoPoints(long n, Photo3 *seq); 877 878 /* --------- Surface Definition ------------------- */ 879 880 /* Compute height of hot spot */ 881 double puffHeight(double h0, double dist); 882 883 /* Compute internal longitude of disk hot spot */ 884 double internal_longitude(); 885 886 /* Map angle in degree into [0, 2*PI) */ 887 double dtor_2pi(double angle); 888 889 /* Define Dimensionless System Parameters */ 890 /* -- Overcontact */ 891 int DefineParamOver(); 892 893 /* Define Dimensionless System Parameters */ 894 int DefineParam(int n); 895 896 /* Define Dimensionless System Parameters for disk */ 897 int DefineParamDisk(); 898 899 /* For power normalization in elliptic orbit */ 900 void ResetNormPower(); 901 902 /* Surface Division - filter dependent */ 903 int DivideSurface(int Comp, int flag_full); 904 905 /* Disk Surface Division - filter dependent */ 906 /* MK 14.03.01 */ 907 int DivideDisk(int Comp); 908 909 /* update radius and potential */ 910 int UpdateGeometry(int n); 911 912 /* --------- Surface Map ------------------- */ 913 914 /* set base path */ 915 int SetMapPath (); 916 917 /* set passband */ 918 void SetMapBand (); 919 920 /* print map */ 921 void PrintSurfaceMap (int thisPhaseStep); 922 923 /* --------- Utilities ------------------- */ 924 925 /* error handler */ 926 void /*@exits@*/ nf_error(char error_text[]); 927 928 /* allocate a matrix */ 929 /*@only@*/ /*@null@*/ /*@out@*/ double **dmatrix(long rows, long columns); 930 931 /* free a matrix */ 932 void free_dmatrix(double **matrix, long rows, long columns); 933 /* convert to lowercase */ 934 void nf_strlwr(char * instring); 935 936 /* clear data file list */ 937 void ClearList(); 938 939 /* add to data file list */ 940 void AddToList(const char * item, int Format); 941 942 /* determine the locale */ 943 const char *guess_the_locale (); 944 945 /* Set user-defined albedo value(s) */ 946 void SetAlbedo (int Comp, const char * str); 947 int SetLimbCoeff(int Comp, const char * str); 948 int isAlbedoSet(int Comp); 949 float GetAlbedo(int Comp); 950 951 /* Compute tickmarks */ 952 void DefineAxis(double *Min, double *Max, double *Inc); 953 954 /**************************************************************************/ 955 /* */ 956 /* GLOBAL VARIABLES */ 957 /* */ 958 /**************************************************************************/ 959 960 #if 0 961 extern char * program_invocation_name; 962 #endif 963 964 /* ----------- Error Messages ---------------------------------------- */ 965 966 extern char *errmsg[23]; 967 968 /* ----------- Flags ------------------------------------------------- */ 969 970 extern FlagsHandle Flags; 971 972 /* ----------- Flux ------------------------------------------------- */ 973 974 /* FluxOut is one phase, filled in LightSynthesize.c::FluxOut */ 975 extern /*@out@*/ PhotoOut *FluxOut; 976 /* FluxOutFull may cover multiple phases if multiple are computed */ 977 extern /*@out@*/ PhotoOut *FluxOutFull; 978 extern int N_FluxOutFull; 979 extern double L90[3][NUM_MAG]; /* Luminosities */ 980 extern double F90[NUM_MAG]; /* Normalization Variable */ 981 extern double D90[3][NUM_MAG]; /* Distance to norm. point*/ 982 extern double P90[NUM_MAG]; /* Normalization Phase */ 983 extern int Umag, Bmag, Vmag; 984 extern int umag, bmag, vmag, ymag; 985 extern int Rmag, Imag; 986 extern int Jmag, Hmag, Kmag; 987 extern int mag1, mag2; 988 extern char Filters[NUM_MAG+2][FILTERNAME_SIZE]; /* the names */ 989 990 /* -----------limb darkening coefficients ------------------------------ */ 991 992 extern double Limb[3][NUM_MAG][4]; /* component, passband, number */ 993 extern double WaveL[3][NUM_MAG]; /* wavelengths */ 994 extern float mono_zero[NUM_MAG]; /* monochromatic wavelengths */ 995 996 /* ----------- Input File(s) -------------------------------------------- */ 997 998 extern Photo3 *Observed[NUM_MAG+2]; 999 extern /*@null@*/ FileType *FileList; 1000 1001 /* ----------- Components ----------------------------------------------- */ 1002 1003 extern const int Primary; /* literal definition */ 1004 extern const int Secondary; /* literal definition */ 1005 extern const int Disk; /* literal definition */ 1006 extern char *Component[]; 1007 #ifdef HAVE_DISK 1008 #define NUM_COMP 3 1009 #else 1010 #define NUM_COMP 2 1011 #endif 1012 extern BinaryComponent Binary[NUM_COMP]; 1013 extern SurfaceElement *Surface[NUM_COMP]; /* stars+disk */ 1014 1015 1016 extern SpotLight *Spot[2]; /* surface spots */ 1017 extern OrbitType Orbit; /* orbit details */ 1018 1019 /* ---------- line profile --------------------------------------------- */ 1020 1021 extern /*@out@*/ float /* **ProData; */ ProData[PHASESTEPS][PROFILE_ARRAY+1]; 1022 1023 /* ---------- misc ---------------------------------------------------- */ 1024 1025 extern char Out_Plot_File[256]; /* output plot file */ 1026 extern int StepsPerPi; /* No of steps per Pi */ 1027 extern double MaximumElements; /* max. no of surface elements*/ 1028 extern int PhaseSteps; /* phasesteps for lightcurve */ 1029 extern int NumPhases; /* how many phases to compute */ 1030 1031 extern double MCMC_param[SDIM]; 1032 extern double MCMC_sdev[SDIM]; 1033 extern int MCMC_num; 1034 1035 /* --- These are Global Variables that are used (mainly) to pass ------- */ 1036 /* --- Arguments to Functions that are called by RootFind/MinFind ------ */ 1037 1038 extern double Mq; /* mass ratio */ 1039 extern double F; /* nonsync rotation */ 1040 extern double RochePot; /* Roche Potential */ 1041 extern double PhaseScaledVolume; /* scaled volume */ 1042 extern double lambda, mu, nu; /* Coordinates */ 1043 1044 1045 /**************************************************************************/ 1046 /* */ 1047 /* INCLUDES */ 1048 /* */ 1049 /**************************************************************************/ 1050 1051 #ifndef HAVE_MKFIFO 1052 #ifdef _WITH_GNUPLOT 1053 #undef _WITH_GNUPLOT 1054 #undef _WITH_PGPLOT 1055 #endif 1056 #endif 1057 1058 #ifndef HAVE_UNISTD_H 1059 1060 #ifdef _WITH_GNUPLOT 1061 #undef _WITH_GNUPLOT 1062 #undef _WITH_PGPLOT 1063 #endif 1064 1065 #endif 1066 1067 /* include limits for parameters */ 1068 1069 #include "LightLimits.h" 1070 1071 1072 /* ---------- MPI ---------------------------------------------------- */ 1073 1074 /* include header for MPI */ 1075 #if defined(_WITH_MPI_COARSE) 1076 void NF_Mpi_distribute (); 1077 void NF_Mpi_collect (); 1078 #endif 1079 1080 #if defined(_WITH_MPI) || defined(_WITH_MPI_COARSE) 1081 1082 #include "mpi.h" 1083 #define NF_MPI_DEBUG 1084 #define NF_MSG_TAG 42 1085 1086 #else 1087 1088 #define MPI_MAX_PROCESSOR_NAME 256 1089 #define MPI_MAX_ERROR_STRING 512 1090 #define MPI_COMM_WORLD 0 1091 /* dummy functions */ 1092 #define MPI_Abort(a, b) 1093 #define MPI_Finalize() 1094 1095 #endif 1096 1097 #define NF_TAG_END 1 1098 #define NF_TAG_MAINLOOP 4 1099 extern int myrank; 1100 extern int numprocs; 1101 extern char processor_name[MPI_MAX_PROCESSOR_NAME]; 1102 extern char error_name[MPI_MAX_ERROR_STRING]; 1103 void NF_Mpi_call(int tag); 1104 void NF_Mpi_endcall(); 1105 void NF_MPI_Abort(); 1106 1107 /* ---------- MPI END ------------------------------------------------- */ 1108 1109 /* include header for PGPLOT emulation */ 1110 1111 #ifdef _WITH_GNUPLOT 1112 #include "LightPGemulate.h" 1113 #endif 1114 1115 /* include header for GUI, also included by gnome.h */ 1116 1117 #ifdef _WITH_GTK 1118 1119 #ifdef HAVE_GNOME 1120 #include <gnome.h> 1121 #else 1122 #include <gtk/gtk.h> 1123 #endif 1124 1125 #include "LightGtk.h" 1126 1127 #ifdef _WITH_OPENGL 1128 /* need to #undef this because /usr/include/jconfig.h 1129 * defines HAVE_STDDEF_H (which is plain stupid IMHO) 1130 */ 1131 #undef HAVE_STDDEF_H 1132 #include "LightGL.h" 1133 #endif 1134 #endif 1135 1136 /* end include header for GUI */ 1137 1138 /* i18n, this stuff is already included by gnome.h */ 1139 1140 #ifndef HAVE_GNOME 1141 1142 #include "gettext.h" 1143 1144 #ifdef ENABLE_NLS 1145 1146 /* #include <libintl.h> */ 1147 #define _(string) gettext (string) 1148 1149 #ifdef gettext_noop 1150 #define N_(string) gettext_noop (string) 1151 #else 1152 #define N_(String) (String) 1153 #endif 1154 1155 #else 1156 1157 #define _(string) string 1158 #define N_(string) string 1159 1160 #endif 1161 1162 /* #ifndef HAVE_GNOME */ 1163 #endif 1164 1165 /* end i18n */ 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176