/* NIGHTFALL Light Curve Synthesis Program */
/* Copyright (C) 1998 Rainer Wichmann */
/* */
/* This program is free software; you can redistribute it */
/* and/or modify */
/* it under the terms of the GNU General Public License as */
/* published by */
/* the Free Software Foundation; either version 2 of the License, or */
/* (at your option) any later version. */
/* */
/* This program is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
/* GNU General Public License for more details. */
/* */
/* You should have received a copy of the GNU General Public License */
/* along with this program; if not, write to the Free Software */
/* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
#include <math.h>
#include <stdio.h>
#include <string.h>
#include <float.h>
#include "Light.h"
#ifdef _WITH_PGPLOT
#ifndef _WITH_GNUPLOT
#include "cpgplot.h"
#endif
#endif
#ifdef _WITH_GTK
#ifdef _WITH_PGPLOT
/* global variables */
static float slicePosition = 0; /* position of slice */
static int slicePlane = 1; /* plane of slice (ON == xy) */
static int sliceState; /* store current plane */
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.1
@short Compute potential in cartesian coordinates
@param (double) x The x position
@param (double) y The y position
@param (double) z The z position
@param (double) q The mass ratio
@param (double) f Asynchronicity
@return (double) The potential
@heading Plotting
****************************************************************************/
double RocheXYSlice(double x, double y, double z, double q, double f)
{
double Eval; /* return value = potential */
double r1, r2; /* radius to point from primary/secondary */
double h1, h2, h3;/* helper variables */
h1 = SQR(x);
h2 = SQR(y);
h3 = SQR(z);
r1 = sqrt(h1 + h2 + h3);
h1 = SQR(x-1);
r2 = sqrt(h1 + h2 + h3);
if (r1 > DBL_EPSILON && r2 > DBL_EPSILON)
{
h1 = SQR(x);
h3 = SQR(f);
Eval = 1.0/r1
+ q * ( 1.0/r2 - x )
+ ( (q+1.0)/2.0 ) * (h1 + h2) * h3;
}
else
{
if (Flags.debug[VERBOSE] == ON || Flags.debug[WARN] == ON)
fprintf(stderr, "**Warning**: Divide by zero, File %s, Line %d \n",
__FILE__, __LINE__);
Eval = 1.0;
}
return(Eval);
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.37
@short Contour plot a slice through the potential
@param (int) GtkEnd What to plot
@return (void)
@heading Plotting
****************************************************************************/
void PlotGtkSlice(int GtkEnd)
{
register unsigned long i, j = 0; /* loop variables */
float RocheXYArray[120][180]; /* matrix for slice */
float x, y; /* position in plane */
float q; /* mass ratio */
float f; /* asynchronicity */
float cont_levels[6]; /* contour levels */
int Ncont; /* number of contour levels */
char sliceTitle[32]; /* plot title */
float transform[6] = {-1.0-0.5*(1.0/60.), /* transformation matrix */
(1.0/60.),
0,
-1.0-0.5*(1.0/60.),
0,
(1.0/60.)};
float *a_ptr = &RocheXYArray[0][0]; /* pointer to slice matrix */
char file_out[256+4]; /* ps output file */
sprintf(file_out, "%s%s", Out_Plot_File, "/CPS");
/* >>>>>>>>>>>>>>>>>>>> open plot device <<<<<<<<<<<<<<<<<<<<<<< */
#ifdef _WITH_GNUPLOT
if (GtkEnd == 3) {
if (cpgopen("/XSERVE") != 1)
{
make_dialog (_(errmsg[2]));
return;
}
}
if (GtkEnd == 2) {
cpgend(); /* end previous interactive plot */
if (cpgopen(file_out) != 1)
{
make_dialog (_(errmsg[1]));
return;
}
}
gnu_start();
#else
if (GtkEnd != 2) {
if(cpgopen("/XSERVE") != 1)
{
make_dialog (_(errmsg[2]));
return;
}
} else {
if(cpgopen(file_out) != 1)
{
make_dialog (_(errmsg[1]));
return;
}
}
#endif
/* >>>>>>>>>>>>>>>>>>>> compute slice <<<<<<<<<<<<<<<<<<<<<<< */
q = Binary[Primary].Mq;
if (Flags.asynchron1 == ON) f = Binary[Primary].Fratio;
else f = 1.0;
for (i = 0; i < 120; ++i) {
for (j = 0; j < 180; ++j) {
y = -1.0 + (i/60.0);
x = -1.0 + (j/60.0);
if (slicePlane == 1)
RocheXYArray[i][j] = -RocheXYSlice(x, y, slicePosition, q, f);
else
RocheXYArray[i][j] = -RocheXYSlice(x, slicePosition, y, q, f);
}
}
cont_levels[0] = -Binary[Primary].RCLag1;
cont_levels[1] = -Binary[Primary].RCLag2;
cont_levels[2] = -Binary[Primary].RCLag3;
cont_levels[3] = -Binary[Primary].RochePot;
RochePot = 0.0;
Mq = Binary[Primary].Mq;
cont_levels[4] = -RochePerpendSecond(Binary[Secondary].Radius);
/* cont_levels[4] = -Binary[Secondary].RochePot; */
Ncont = 5;
/* >>>>>>>>>>>>>>>>>>>> plot slice <<<<<<<<<<<<<<<<<<<<<<< */
cpgscf(2); cpgslw(1); cpgsch(1.2);
cpgsvp(0.15,0.9,0.1,1.0);
cpgwnad(-1., 2., -1.0, 1.0);
if (slicePlane == 1) {
/* slice is parallel to xy plane */
sprintf(sliceTitle, _("Roche Potential at Z = %4.2f"), slicePosition);
cpglab(_("X"), _("Y"), sliceTitle);
} else {
/* slice is parallel to xz plane */
sprintf(sliceTitle, _("Roche Potential at Y = %4.2f"), slicePosition);
cpglab(_("X"), _("Z"), sliceTitle);
}
#ifdef _WITH_GNUPLOT
cpgbox("BCNST", 0.5, 0, "BCNST", 0.5, 0);
#else
cpgbox("BCNST", 0,0, "BCNST", 0,0);
#endif
cpgcont(a_ptr, 180, 120, 1, 180, 1, 120,
cont_levels, (-Ncont), transform);
/* >>>>>>>>>>>>>>>>>>>> end plot <<<<<<<<<<<<<<<<<<<<<<< */
if (GtkEnd == 2) my_cpgend();
else cpgend();
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Make hardcopy of slice (ps file)
@param (GtkWidget) *widget Discarded
@param (gpointer) *data Discarded
@return (void)
@heading Plotting
****************************************************************************/
void hard_slice (GtkWidget *widget, gpointer *data)
{
PlotGtkSlice (2);
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Destroy the Slicer
@param (GtkWidget) *widget Discarded
@param (gpointer) *data Widget to destroy
@return (void)
@heading Plotting
****************************************************************************/
void delete_slice (GtkWidget *widget, gpointer *data)
{
gtk_widget_destroy (GTK_WIDGET (data) );
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Adjust slice position
@param (GtkAdjustment) *adj The new position
@param (gpointer) *data Discarded
@return (void)
@heading Plotting
****************************************************************************/
void slice_change (GtkAdjustment *adj, GtkWidget *data)
{
slicePosition = 0.01 * adj->value;
PlotGtkSlice (ON);
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Change slice plane
@param (GtkWidget) *widget Discarded
@param (gpointer) *data Discarded
@return (void)
@heading Plotting
****************************************************************************/
void slice_plane (GtkWidget *widget, gpointer *data)
{
/* 0 = xz, 1 = xy */
sscanf ((char *) data, "%d", &slicePlane);
if (slicePlane != sliceState) {
PlotGtkSlice (ON);
sliceState = slicePlane;
}
return;
}
/****************************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short Interactive display of slices through Roche potential
@param (void)
@return (void)
@heading Plotting
****************************************************************************/
void MakeSliceBox()
{
GtkWidget *label;
GtkWidget *button;
GtkWidget *vbox;
GtkWidget *hbox;
GtkWidget *frame;
GtkWidget *separator;
GtkWidget *scale1;
GtkObject *adj1;
GtkWidget *slice_window;
GSList *group;
GtkTooltips *tooltips;
/* >>>>>>>>>>>>>>>>>>>>>>> MAKE WIDGET <<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
tooltips = gtk_tooltips_new ();
gtk_tooltips_set_delay (tooltips, 1200);
slice_window = gtk_window_new (GTK_WINDOW_TOPLEVEL);
/* gtk_widget_set_usize (slice_window, 200, 150); */
gtk_window_set_position(GTK_WINDOW(slice_window), GTK_WIN_POS_MOUSE);
#ifdef USING_GTK2
gtk_window_set_resizable(GTK_WINDOW(slice_window), TRUE);
#else
gtk_window_set_policy (GTK_WINDOW(slice_window), FALSE, FALSE, TRUE);
#endif
/* gtk_window_set_policy (GTK_WINDOW(slice_window), TRUE, TRUE, FALSE); */
gtk_signal_connect (GTK_OBJECT (slice_window), "destroy",
(GtkSignalFunc) delete_slice,
GTK_WIDGET (slice_window));
gtk_window_set_title (GTK_WINDOW (slice_window), _("Roche Slicer") );
gtk_container_set_border_width (GTK_CONTAINER (slice_window), 0);
frame = gtk_frame_new (_("Viewing Options") );
gtk_container_add (GTK_CONTAINER (slice_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);
adj1 = gtk_adjustment_new (0.0, -100.0, 100.0, 0.1, 1.0, 1.0);
gtk_signal_connect (GTK_OBJECT (adj1), "value_changed",
GTK_SIGNAL_FUNC (slice_change), 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);
scale1 = gtk_hscale_new(GTK_ADJUSTMENT (adj1));
#ifdef _WITH_GNUPLOT
/* with GTK_UPDATE_CONTINUOUS, I get the error: */
/* ** ERROR **: sigpipe caught */
/* Broken pipe */
/* is contour plotting in gnuplot too slow ?? or the fprintf to fifo ?? */
gtk_range_set_update_policy (GTK_RANGE (scale1), GTK_UPDATE_DISCONTINUOUS);
#else
gtk_range_set_update_policy (GTK_RANGE (scale1), GTK_UPDATE_CONTINUOUS);
#endif
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 (_("Slice Position") );
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);
button = gtk_radio_button_new_with_label (NULL, _("X-Z Plane") );
gtk_signal_connect (GTK_OBJECT (button), "clicked",
(GtkSignalFunc) slice_plane,
(gpointer) "0");
gtk_box_pack_start (GTK_BOX (hbox), button, TRUE, TRUE, 0);
gtk_widget_show (button);
group = gtk_radio_button_group (GTK_RADIO_BUTTON (button));
button = gtk_radio_button_new_with_label(group, _("X-Y Plane") );
gtk_toggle_button_set_active (GTK_TOGGLE_BUTTON (button), TRUE);
gtk_signal_connect (GTK_OBJECT (button), "clicked",
(GtkSignalFunc) slice_plane,
(gpointer) "1");
gtk_box_pack_start (GTK_BOX (hbox), button, TRUE, TRUE, 0);
gtk_widget_show (button);
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_slice,
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_slice,
GTK_WIDGET (slice_window));
gtk_tooltips_set_tip (tooltips, button,
_("Close Roche Slicer"),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 (slice_window);
/* >>>>>>>>>>>>>>>>>>>>>>> initial plot <<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
slicePlane = 1; sliceState = 1;
PlotGtkSlice (3);
return;
}
/* with pgplot */
#endif
/* with gtk */
#endif
syntax highlighted by Code2HTML, v. 0.9.1