/* NIGHTFALL Light Curve Synthesis Program */ /* Copyright (C) 1998 Rainer Wichmann */ /* */ /* This program is free software; you can redistribute it */ /* and/or modify */ /* it under the terms of the GNU General Public License as */ /* published by */ /* the Free Software Foundation; either version 2 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program; if not, write to the Free Software */ /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */ #include #include #include #include #include "Light.h" #ifdef _WITH_PGPLOT #ifndef _WITH_GNUPLOT #include "cpgplot.h" #endif #endif #ifdef _WITH_GTK #ifdef _WITH_PGPLOT /************************************************************* Variables for temporary storage of data **************************************************************/ static double vis_keepInclination; /* store current inclination */ static double vis_keepPhase; /* store current phase */ static int visBlock = ON; /* block input */ /****************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Simplified version of MainLoop() to compute visibility @param (int) JJ The phase index @return (int) The exit status @heading Plotting *******************************************************************/ int MainSingle(int JJ) { int PhaseIndex = JJ; /* phase index */ int ErrCode; /* exit status subroutines */ register unsigned long i; /* loop variable */ SurfaceElement *SurfPtrP; /* pointer to surface primary */ SurfaceElement *SurfPtrS; /* pointer to surface secondary */ #ifdef HAVE_DISK SurfaceElement *SurfPtrD; /* pointer to surface disk */ #endif SurfPtrP = Surface[Primary]; SurfPtrS = Surface[Secondary]; #ifdef HAVE_DISK SurfPtrD = Surface[Disk]; #endif for (i = 0; i < MAXELEMENTS; ++i) { SurfPtrP->SpotNum = 0; SurfPtrS->SpotNum = 0; SurfPtrP->SumDimFactor = 0; SurfPtrS->SumDimFactor = 0; ++SurfPtrP; ++SurfPtrS; } /* >>>>>>>>>>>>> LightGeometry <<<<<<<<<<<<<<< */ Flags.first_pass = ON; check_data(); if (Flags.elliptic == ON) FindPeriastron(); if (Flags.fill == OFF) { ErrCode = DefineParam(Primary); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) printf(_("\n -- Primary defined --\n")); DefineParam(Secondary); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) printf(_("\n -- Secondary defined --\n")); #ifdef HAVE_DISK if (Flags.disk == ON) { DefineParamDisk(); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) printf(_("\n -- Disk defined --\n")); } #endif } else { ErrCode = DefineParamOver(); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) printf(_("\n -- Overcontact Geometry defined --\n")); } /* >>>>>>>>>>>>> LightDivide <<<<<<<<<<<<<<<< */ ErrCode = DivideSurface(Primary); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) { printf(_("\n Primary divided\n") ); } ErrCode = DivideSurface(Secondary); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) printf(_("\n Secondary divided\n") ); #ifdef HAVE_DISK if (Flags.disk == ON) { ErrCode = DivideDisk(Disk); if (ErrCode == 8) return(8); if (Flags.debug[VERBOSE] == ON) printf(_("\n Disk divided\n") ); } #endif /* >>>>>>>>>>>>> test configuration <<<<<<<<<<<<<<< */ if(LightSetupTests() > 0) { if (Flags.interactive == OFF) { nf_error(_(" --> Program Aborted")); } else { WARNING(_(" --> Program Aborted")); return(8); } } Flags.first_pass = OFF; /* ------------- LightElliptic --------------------- */ Orbit.PhaseIndex = PhaseIndex; if (Flags.elliptic == ON) { /* update (nearly) everything */ SolveKepler(PhaseIndex); ErrCode = UpdateGeometry(Primary); if (ErrCode == 8) return(8); ErrCode = UpdateGeometry(Secondary); if (ErrCode == 8) return(8); /* ------- IT IS IMPERATIVE TO DO SECONDARY SECOND !!! ---- */ ErrCode = DivideSurface(Primary); if (ErrCode == 8) return(8); ErrCode = DivideSurface(Secondary); if (ErrCode == 8) return(8); ErrCode = LightSetupTests(); if (ErrCode > 0) return(8); } /* ----------- move spots if asynchroneous rotation ----------- */ if ( (Flags.Spots1 > OFF) && (Flags.asynchron1 == ON && Flags.elliptic == OFF) ) { MakeSpots(Primary, Orbit.PhaseIndex); } if ( (Flags.Spots2 > OFF) && ( Flags.asynchron2 == ON && Flags.elliptic == OFF) ) { MakeSpots(Secondary, Orbit.PhaseIndex); } /* >>>>>>>>>>>>> LightSynthesize <<<<<<<<<<<<<<< */ if (Flags.fill == ON) LightCopyThroat(); LightCurve(PhaseIndex); #ifdef HAVE_DISK if (Flags.disk == ON) LightCurveDisk(PhaseIndex); #endif if (Flags.fill == ON) LightCopyThroatBack(); if (Flags.elliptic == ON) PhaseShift(); return(0); } /****************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Plot a view of the binary system @param (float) GetPhase The phase @param (float) GetInclination The inclination @param (int) GtkEnd A Flag @return (void) @heading Plotting *******************************************************************/ void PlotGtkGeometry(float GetPhase, float GetInclination, int GtkEnd) { float GtkInclination = Orbit.Inclination; /* inclination */ float GtkPhase; /* phase */ register unsigned long i, JJ = 0; /* loop variables */ long CountElem; /* # of elements */ long Comp; /* stellar component */ long MaxComp = 2; /* # of components */ float *Xplot; float *Yplot; double sinAlpha; /* sin inclination */ double cosAlpha; /* cos inclination */ double sinPhi; /* sin phase */ double cosPhi; /* cos phase */ double Trans[4][4]; /* transformation matrix */ double t10, t11, t12, t13; /* transformation matrix */ double t20, t21, t22, t23; /* transformation matrix */ double C; /* sign (for primary/secondary) */ SurfaceElement *SurfPtr; /* pointer to surface */ int ErrCode; /* exit status */ char visTitle[32]; /* plot title */ /* ----------- initialize ----------------------------------- */ double storeInclination = Orbit.Inclination, storePhase = Orbit.Phase; char file_out[256+4]; /* ps output file */ sprintf(file_out, "%s%s", Out_Plot_File, "/CPS"); Orbit.Inclination = vis_keepInclination; JJ = vis_keepPhase; /* ----------- allocate memory and check -------------------- */ Xplot = malloc(2*MAXELEMENTS*sizeof(float)); Yplot = malloc(2*MAXELEMENTS*sizeof(float)); if (Xplot == NULL || Yplot == NULL) { if (Flags.interactive == ON) { if (Xplot != NULL) free(Xplot); if (Yplot != NULL) free(Yplot); make_dialog(_(errmsg[0])); return; } else nf_error(_(errmsg[0])); } /* >>>>>>>>>>>>>>>>>> BEGIN GEOMETRY <<<<<<<<<<<<<<<<<< */ /* ------------- set variables ----------------------- */ if (GtkEnd == 0 || GtkEnd == -2) { JJ = PhaseSteps * GetPhase; if (JJ < 0) JJ = 0.; if (JJ >= PhaseSteps) JJ = PhaseSteps-1; vis_keepPhase = JJ; Orbit.Inclination = vis_keepInclination; } if (GtkEnd == 0 || GtkEnd == -1) { Orbit.Inclination = GetInclination; vis_keepInclination = Orbit.Inclination; } /* --------------- validate input ---------------------- */ check_data(); /* --------------- visibility ---------------------- */ ErrCode = MainSingle(JJ); if (ErrCode != 0) { free(Xplot); free(Yplot); make_dialog (_(errmsg[12])); Orbit.Phase = storePhase; Orbit.Inclination = storeInclination; return; } /* >>>>>>>>>>>>>>>>>>>>>>> END GEOMETRY <<<<<<<<<<<<<<<<<<<<< */ GtkInclination = Orbit.Inclination; if (Flags.elliptic == OFF) { GtkPhase = Orbit.Phase; if ( -0.5 + (GtkPhase/(2.0*PI)) >= 0.0) sprintf(visTitle, _("Phase %5.2f"), -0.5 + (GtkPhase/(2.0*PI)) ); else sprintf(visTitle, _("Phase %5.2f"), -0.5 + (GtkPhase/(2.0*PI)) + 1.0); } else { GtkPhase = Orbit.Phase; if (-0.5 + ((Orbit.M + PI + Orbit.OmegaZero)/(2.0*PI)) >= 0.0) sprintf(visTitle, _("Phase %5.2f"), -0.5 + ((Orbit.M + PI + Orbit.OmegaZero)/(2.0*PI)) ); else sprintf(visTitle, _("Phase %5.2f"), -0.5 + ((Orbit.M + PI + Orbit.OmegaZero)/(2.0*PI)) + 1.0); } #ifdef _WITH_GNUPLOT if (GtkEnd == 0 || GtkEnd == 3) { if (cpgopen("/XSERVE") != 1) { make_dialog (_(errmsg[2])); free(Xplot); free(Yplot); Orbit.Phase = storePhase; Orbit.Inclination = storeInclination; return; } } if (GtkEnd == 2) { if (cpgopen(file_out) != 1) { make_dialog (_(errmsg[1])); free(Xplot); free(Yplot); Orbit.Phase = storePhase; Orbit.Inclination = storeInclination; return; } } gnu_start(); #else if (GtkEnd != 2) { if(cpgopen("/XSERVE") != 1) { make_dialog (_(errmsg[2])); free(Xplot); free(Yplot); Orbit.Phase = storePhase; Orbit.Inclination = storeInclination; return; } } else { if(cpgopen(file_out) != 1) { make_dialog (_(errmsg[1])); free(Xplot); free(Yplot); Orbit.Phase = storePhase; Orbit.Inclination = storeInclination; return; } } #endif sinAlpha = sin(GtkInclination - 0.5*PI); cosAlpha = cos(GtkInclination - 0.5*PI); /* plot box */ cpgscf(2); cpgslw(1); cpgsch(0.9); cpgsvp(0.15 ,0.85 ,0.15 , 0.85); cpgwnad(-1.2, 1.2, -0.98, 0.98); #ifdef _WITH_GNUPLOT cpgbox("BCNST", 0,0, "BCNST", 0,0); #endif cpgsch(2.0); cpgmtxt("T", 1, 0.8, 0.5, visTitle); cpgsch(0.9); /* reset size */ cpgsch(0.1); cpgslw(1); /* loop over components */ CountElem = 0; #ifdef HAVE_DISK if (Flags.disk == ON) MaxComp = 3; #endif for (Comp = 0; Comp < MaxComp; ++Comp) { /* now we need the phase */ if (Comp == Primary) { sinPhi = sin(-GtkPhase); cosPhi = cos(-GtkPhase); C = 1.0; } else { C = -1.0; /* mirror y -> (-y) for secondary */ sinPhi = sin(-GtkPhase + PI); cosPhi = cos(-GtkPhase + PI); } /* the transformation matrix */ /* first column, then row */ Trans[0][0] = cosPhi * cosAlpha; Trans[1][0] = sinPhi; Trans[2][0] = cosPhi * sinAlpha; Trans[3][0] = 0.0; Trans[0][1] = -sinPhi * cosAlpha; Trans[1][1] = cosPhi; Trans[2][1] = -sinPhi * sinAlpha; Trans[3][1] = 0.0; Trans[0][2] = -sinAlpha; Trans[1][2] = 0.0; Trans[2][2] = cosAlpha; Trans[3][2] = 0.0; Trans[0][3] = Binary[Comp].RLag1 * (1.0 - cosPhi * cosAlpha); Trans[1][3] = -Binary[Comp].RLag1 * sinPhi; Trans[2][3] = -Binary[Comp].RLag1 * cosPhi * sinAlpha; Trans[3][3] = 1.0; t10 = Trans[1][0]; t11 = Trans[1][1]; t12 = Trans[1][2]; t13 = Trans[1][3]; t20 = Trans[2][0]; t21 = Trans[2][1]; t22 = Trans[2][2]; t23 = Trans[2][3]; /* do the transformation */ SurfPtr = Surface[Comp]; for (i = 0; i < Binary[Comp].NumElem; ++i) { if (SurfPtr->EclipseFlag == 0 && SurfPtr->SpotNum == 0) { /* Observer looks from positive X-Axis */ /* this is the negative Y-Axis */ Xplot[CountElem] = - (SurfPtr->lambda * t10 + SurfPtr->mu * t11 * C + SurfPtr->nu * t12 + 1.0 * t13); /* and this the Z-Axis */ Yplot[CountElem] = (SurfPtr->lambda * t20 + SurfPtr->mu * t21 * C + SurfPtr->nu * t22 + 1.0 * t23); /* increase counter for elements in Xplot/Yplot */ ++CountElem; } ++SurfPtr; } } /* plot the data */ #ifdef _WITH_GNUPLOT cpgsvp(0.05,0.9,0.05,0.95); cpgwnad(-1.2, 1.2, -0.9, 0.9); cpgmtxt("T", 1, 0.8, 0.5, visTitle); #endif cpgpt(CountElem, Xplot, Yplot, 17); /* plot finished */ if (GtkEnd == 2) my_cpgend(); else cpgend(); Orbit.Phase = storePhase; Orbit.Inclination = storeInclination; free(Xplot); free(Yplot); return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Make hardcopy of StarView (ps file) @param (GtkWidget) *widget Discarded @param (gpointer) *data Discarded @return (void) @heading Plotting ****************************************************************************/ void hard_vis (GtkWidget *widget, gpointer data) { if (visBlock == OFF) PlotGtkGeometry (0.0, 0.0, 2); return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Destroy StarView @param (GtkWidget) *widget Discarded @param (gpointer) *data Widget to destroy @return (void) @heading Plotting ****************************************************************************/ void delete_vis (GtkWidget *widget, gpointer data) { if (visBlock == OFF) PlotGtkGeometry(0.0, 0.0, 1); visBlock = ON; gtk_widget_destroy (GTK_WIDGET (data) ); Binary[Primary].RocheFill = Binary[Primary].RocheStore; Binary[Secondary].RocheFill = Binary[Secondary].RocheStore; #ifdef HAVE_DISK Binary[Disk].RocheFill = Binary[Disk].RocheStore; #endif return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Change the inclination @param (GtkAdjustment) *adj The inclination @param (gpointer) *data Discarded @return (void) @heading Plotting ****************************************************************************/ void vis_rotate_incl (GtkAdjustment *adj, GtkWidget *data) { if (visBlock == OFF) PlotGtkGeometry(0.0, DTOR*(adj->value), -1); Binary[Primary].RocheFill = Binary[Primary].RocheStore; Binary[Secondary].RocheFill = Binary[Secondary].RocheStore; #ifdef HAVE_DISK Binary[Disk].RocheFill = Binary[Disk].RocheStore; #endif return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Change the phase @param (GtkAdjustment) *adj The phase @param (gpointer) *data Discarded @return (void) @heading Plotting ****************************************************************************/ void vis_rotate_phase (GtkAdjustment *adj, GtkWidget *data) { if (visBlock == OFF) PlotGtkGeometry((adj->value / 360.), 0.0, -2); Binary[Primary].RocheFill = Binary[Primary].RocheStore; Binary[Secondary].RocheFill = Binary[Secondary].RocheStore; #ifdef HAVE_DISK Binary[Disk].RocheFill = Binary[Disk].RocheStore; #endif return; } /**************************************************************************** @package nightfall @author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de) @version 1.0 @short Interactive visualization of stars @param (void) @return (void) @heading Plotting ****************************************************************************/ void MakeVbox() { GtkWidget *label; GtkWidget *button; GtkWidget *vbox; GtkWidget *hbox; GtkWidget *separator; GtkWidget *scale1; GtkWidget *scale2; GtkWidget *frame; GtkObject *adj1, *adj2; GtkWidget *vis_window; GtkTooltips *tooltips; /* ----------- initialize ---------------------- */ tooltips = gtk_tooltips_new (); gtk_tooltips_set_delay (tooltips, 1200); Binary[Primary].RocheStore = Binary[Primary].RocheFill; Binary[Secondary].RocheStore = Binary[Secondary].RocheFill; #ifdef HAVE_DISK Binary[Disk].RocheStore = Binary[Disk].RocheFill; #endif /* >>>>>>>>>>>>>>>>>>>> MAKE WIDGET <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ vis_window = gtk_window_new (GTK_WINDOW_TOPLEVEL); /* gtk_widget_set_usize (vis_window, 200, 170); */ gtk_window_set_position(GTK_WINDOW(vis_window), GTK_WIN_POS_MOUSE); #ifdef USING_GTK2 gtk_window_set_resizable(GTK_WINDOW(vis_window), TRUE); #else gtk_window_set_policy (GTK_WINDOW(vis_window), FALSE, FALSE, TRUE); #endif /* gtk_window_set_policy (GTK_WINDOW(vis_window), TRUE, TRUE, FALSE); */ gtk_signal_connect (GTK_OBJECT (vis_window), "destroy", (GtkSignalFunc) delete_vis, GTK_WIDGET (vis_window)); gtk_window_set_title (GTK_WINDOW (vis_window), _("StarView") ); gtk_container_set_border_width (GTK_CONTAINER (vis_window), 2); frame = gtk_frame_new (_("Viewing Options") ); gtk_container_add (GTK_CONTAINER (vis_window), frame); gtk_container_set_border_width(GTK_CONTAINER(frame), 2); gtk_widget_show (frame); vbox = gtk_vbox_new (FALSE, 0); gtk_container_add (GTK_CONTAINER (frame), vbox); gtk_widget_show (vbox); adj2 = gtk_adjustment_new (RTOD*Orbit.Inclination, 0.0, 90.0, 1.0, 1.0, 1.0); gtk_signal_connect (GTK_OBJECT (adj2), "value_changed", GTK_SIGNAL_FUNC (vis_rotate_incl), NULL); adj1 = gtk_adjustment_new (0.0, 0.0, 360.0, 1.0, 1.0, 1.0); gtk_signal_connect (GTK_OBJECT (adj1), "value_changed", GTK_SIGNAL_FUNC (vis_rotate_phase), NULL); hbox = gtk_hbox_new (TRUE, 0); gtk_container_set_border_width (GTK_CONTAINER (hbox), 4); gtk_box_pack_start (GTK_BOX (vbox), hbox, TRUE, TRUE, 0); gtk_widget_show (hbox); scale2 = gtk_hscale_new(GTK_ADJUSTMENT (adj2)); gtk_range_set_update_policy (GTK_RANGE (scale2), GTK_UPDATE_DELAYED); gtk_box_pack_start (GTK_BOX (hbox), scale2, TRUE, TRUE, 0); gtk_widget_show (scale2); hbox = gtk_hbox_new (TRUE, 0); gtk_container_set_border_width (GTK_CONTAINER (hbox), 4); gtk_box_pack_start (GTK_BOX (vbox), hbox, TRUE, TRUE, 0); gtk_widget_show (hbox); label = gtk_label_new (_("Inclination") ); gtk_box_pack_start (GTK_BOX (hbox), label, TRUE, TRUE, 0); gtk_widget_show (label); hbox = gtk_hbox_new (TRUE, 0); gtk_container_set_border_width (GTK_CONTAINER (hbox), 4); gtk_box_pack_start (GTK_BOX (vbox), hbox, TRUE, TRUE, 0); gtk_widget_show (hbox); scale1 = gtk_hscale_new(GTK_ADJUSTMENT (adj1)); gtk_range_set_update_policy (GTK_RANGE (scale1), GTK_UPDATE_DELAYED); gtk_box_pack_start (GTK_BOX (hbox), scale1, TRUE, TRUE, 0); gtk_widget_show (scale1); hbox = gtk_hbox_new (TRUE, 0); gtk_container_set_border_width (GTK_CONTAINER (hbox), 4); gtk_box_pack_start (GTK_BOX (vbox), hbox, TRUE, TRUE, 0); gtk_widget_show (hbox); label = gtk_label_new (_("Phase Angle") ); gtk_box_pack_start (GTK_BOX (hbox), label, TRUE, TRUE, 0); gtk_widget_show (label); separator = gtk_hseparator_new (); gtk_box_pack_start (GTK_BOX (vbox), separator, FALSE, TRUE, 0); gtk_widget_show (separator); hbox = gtk_hbox_new (FALSE, 0); gtk_box_pack_start (GTK_BOX (vbox), hbox, FALSE, TRUE, 0); gtk_widget_show (hbox); #ifdef USING_GTK2 button = gtk_button_new_with_label (_("Postscript") ); gtk_button_set_image (GTK_BUTTON(button), gtk_image_new_from_stock (GTK_STOCK_PRINT, GTK_ICON_SIZE_SMALL_TOOLBAR) ); #else button = gtk_button_new_with_label (_("Postscript") ); #endif gtk_signal_connect (GTK_OBJECT (button), "clicked", (GtkSignalFunc) hard_vis, NULL); gtk_tooltips_set_tip (tooltips, button, _("Make postscript hardcopy from current plot"), NULL); gtk_box_pack_start (GTK_BOX (hbox), button, TRUE, TRUE, 0); (GTK_WIDGET_FLAGS (button) |= (GTK_CAN_DEFAULT)); gtk_widget_grab_default (button); gtk_widget_show (button); #ifdef USING_GTK2 button = gtk_button_new_from_stock(GTK_STOCK_OK); #else button = gtk_button_new_with_label (_("Close") ); #endif gtk_signal_connect (GTK_OBJECT (button), "clicked", (GtkSignalFunc) delete_vis, GTK_WIDGET (vis_window)); gtk_tooltips_set_tip (tooltips, button, _("Close StarView") ,NULL); gtk_box_pack_start (GTK_BOX (hbox), button, TRUE, TRUE, 0); (GTK_WIDGET_FLAGS (button) |= (GTK_CAN_DEFAULT)); gtk_widget_grab_default (button); gtk_widget_show (button); gtk_widget_show (vis_window); /* ----------- initial plot ---------------------- */ vis_keepPhase = 0; vis_keepInclination = Orbit.Inclination; visBlock = OFF; PlotGtkGeometry(0, Orbit.Inclination, 0); return; } /* with pgplot */ #endif /* with gtk */ #endif