/* 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"
#include "LightSimplex.h"
#ifdef _WITH_PGPLOT
#ifndef _WITH_GNUPLOT
#include "cpgplot.h"
#endif
#endif
#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif
/*********************************************************************
@package nightfall
@author Rainer Wichmann (rwichman@lsw.uni-heidelberg.de)
@version 1.0
@short 2D chi-square map
@param (void)
@return (int) status The exit status
@heading Chi Square Map
**********************************************************************/
int ChiMap()
{
FILE *outfile; /* the output file */
long dof; /* degrees of freedom */
int ErrCode; /* return status of subroutines */
double Merit = 0; /* chi square */
int Count, j, i, k; /* loop variables */
int Test=0; /* test variable */
double Start[2]; /* start values */
double Step[2]; /* step values */
double Scan[SDIM]; /* some Vertex */
double X[SDIM] = { 0.0 }; /* some Vertex */
static float Y[_CHI_SCANS_][ _CHI_SCANS_]; /* map array */
static float P[_CHI_SCANS_][ _CHI_SCANS_]; /* plot array */
int VertexTest; /* test whether within limits */
float Min=0.0, Max=0.0; /* plot scaling */
int fir_p = 0, sec_p = 0; /* fit parameters */
char message[256]; /* status message */
#ifdef _WITH_PGPLOT
float MedTest; /* plot scaling */
float transform[6] = {0, 1, 0, 0, 0, 1};
/* transformation matrix */
float cont_levels[6]; /* contour levels */
int Ncont; /* # of contour levels */
float *a_ptr = &Y[0][0]; /* pointer to plot array */
char file_out[256+4]; /* ps output file */
char cont_labels[9][9]; /* contour labels */
sprintf(file_out, "%s%s", Out_Plot_File, "/CPS");
#endif
/* >>>>>>>>>>>> Initialize <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
dof = 0;
for (j = 0; j < NUM_MAG+2; ++j) {
dof = dof + Flags.Passbands[j];
}
if (dof <= 0) {
#ifdef _WITH_GTK
if ( Flags.interactive == ON) {
make_dialog(_(errmsg[18]));
} else {
WARNING("\n");
WARNING(_(" Dear User,\n"));
WARNING(_(" We have a minor problem: "));
WARNING(_(" There are no data, thus we abort"));
WARNING(_(" the mapping procedure"));
}
return(8);
#else
WARNING("\n");
WARNING(_(" Dear User,\n"));
WARNING(_(" We have a minor problem: "));
WARNING(_(" There are no data, thus we abort"));
WARNING(_(" the mapping procedure"));
return(8);
#endif
}
SimplexInitRanges();
for (i = 0; i < SDIM; ++i) {
Scan[i] = (*Param[i].par);
X[i] = (double) 0.0;
}
for (i = 0; i < _CHI_SCANS_; ++i) {
for (j = 0; j < _CHI_SCANS_; ++j) Y[i][j] = 0.0;
}
/* count parameters to do map for */
Count = 0;
for (j = 0; j < SDIM; ++j) {
if (Flags.simplex[j] == ON) {
if (Count == 0) fir_p = j;
else if (Count == 1) sec_p = j;
else if (Count >= 2) {
Flags.simplex[j] = OFF;
WARNING(_("Too much Parameters for ChiSquare Map, some ignored"));
}
++Count;
}
}
/* error: too few parameters */
if (Count < 2) {
#ifdef _WITH_GTK
if (Flags.interactive == OFF) {
nf_error(_(errmsg[6]));
} else {
make_dialog(_(errmsg[6]));
return(8);
}
#else
nf_error(_(errmsg[6]));
#endif
}
#ifdef _WITH_GTK
if (Flags.interactive == ON)
my_appbar_push(_("Starting Chi Square mapping"));
#endif
Start[0] = 0.0; Start[1] = 0.0;
Step[0] = 1.0; Step[1] = 1.0;
/* initialize start values */
/* this will find the first active parameter */
for (i = SDIM-1; i >= 0; --i) {
if (Flags.simplex[i] == ON) {
Start[0] = Scan[i]; Test = i;
if (i == 1) Step[0] = ( DTOR * Flags.Step[0]);
else if (i == 27) Step[0] = ( 1.989E30 * Flags.Step[0]);
else if (i == 28) Step[0] = ( 6.960E8 * Flags.Step[0]);
else Step[0] = Flags.Step[0];
}
}
/* and this the second */
for (i = SDIM-1; i >= 0; --i) {
if (Flags.simplex[i] == ON && Test != i) {
Start[1] = Scan[i];
if (i == 1) Step[1] = ( DTOR * Flags.Step[1]);
else if (i == 27) Step[1] = ( 1.989E30 * Flags.Step[1]);
else if (i == 28) Step[1] = ( 6.960E8 * Flags.Step[1]);
else Step[1] = Flags.Step[1];
}
}
/* >>>>>>>>>>>>>>>>> CREATE MAP <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
for (i = 0; i < _CHI_SCANS_; ++i) {
for (j = 0; j < _CHI_SCANS_; ++j) {
Count = 0;
for (k = 0; k < SDIM; ++k) {
X[k] = Scan[k];
if (Flags.simplex[k] == ON) {
if (Count == 0) {
X[k] = X[k] + i * Step[0];
}
if (Count == 1) {
X[k] = X[k] + j * Step[1];
}
/* fprintf(stderr, "%8.4g %8.4g\n", X[0], X[1]); */
++Count;
}
}
VertexTest = SimplexCheckVertex(X);
if (VertexTest > 0) WARNING(_("Out of Range, Map will be corrupted"));
SimplexSetVariables(X);
ErrCode = MainLoop(&Merit);
Y[i][j] = (float)(Merit);
P[j][i] = Y[i][j];
if (i == 0 && j == 0) {
Min=Y[i][j];
Max=Y[i][j];
} else {
Min=MIN(Min, Y[i][j]);
Max=MAX(Max,Y[i][j]);
}
}
sprintf(message, _("Column %5d of %5d finished\n"), i+1, _CHI_SCANS_);
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON)
printf(message);
#ifdef _WITH_GTK
if (Flags.interactive == ON) {
sprintf(message, _("Column %5d of %5d finished"), i+1, _CHI_SCANS_);
my_appbar_push(message);
sleep (1);
}
#endif
}
/* >>>>>>>>>>>>>>>>> transform start/step <<<<<<<<<<<<<<<<<<<< */
if (fir_p == 1) {
Step[0] = Flags.Step[0]; Start[0] = Start[0] * RTOD;
}
else if (fir_p == 27) {
Step[0] = Flags.Step[0]; Start[0] = Start[0] / 1.989E30;
}
else if (fir_p == 28) {
Step[0] = Flags.Step[0]; Start[0] = Start[0] / 6.960E8;
}
if (sec_p == 1) {
Step[1] = Flags.Step[1]; Start[1] = Start[1] * RTOD;
}
else if (sec_p == 27) {
Step[1] = Flags.Step[1]; Start[1] = Start[1] / 1.989E30;
}
else if (sec_p == 28) {
Step[1] = Flags.Step[1]; Start[1] = Start[1] / 6.960E8;
}
/* >>>>>>>>>>>>>>>>> PRINT OUT MAP <<<<<<<<<<<<<<<<<<<<<<<<<<< */
if (Flags.debug[BUSY] == ON || Flags.debug[VERBOSE] == ON) {
printf(_(" Chi Square Map:\n"));
/* i = first param, j = second param */
printf(" ");
for (i = 0; i < _CHI_SCANS_; ++i) printf(" %8.4g", Start[0] + i * Step[0]);
printf("\n");
for (j = 0; j < _CHI_SCANS_; ++j) {
printf("%8.4g", Start[1] + j * Step[1]);
for (i = 0; i < _CHI_SCANS_; ++i) printf(" %8.2f", Y[i][j]);
printf("\n");
}
printf(_(" -- end map\n"));
}
outfile = fopen(OUT_CHIMAP, "w");
if (outfile == NULL)
#ifdef _WITH_GTK
{
if ( Flags.interactive == ON) {
make_dialog(_(errmsg[3]));
return (8);
} else nf_error(_(errmsg[3]));
}
#else
nf_error(_(errmsg[3]));
#endif
/* ----------------- write header ------------------------------ */
OutfileHeader(outfile);
fprintf(outfile, _(" Chi Square Map:\n"));
/* i = first param, j = second param */
fprintf(outfile, " ");
for (i = 0; i < _CHI_SCANS_; ++i)
fprintf(outfile, " %8.4g", Start[0] + i * Step[0]);
fprintf(outfile, "\n");
for (j = 0; j < _CHI_SCANS_; ++j) {
fprintf(outfile, "%8.4g", Start[1] + j * Step[1]);
for (i = 0; i < _CHI_SCANS_; ++i) fprintf(outfile, " %8.2f", Y[i][j]);
fprintf(outfile, "\n");
}
fprintf(outfile, _(" -- end map\n"));
/* >>>>>>>>>>>>>>> PLOT MAP <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */
#ifdef _WITH_PGPLOT
transform[0] = Start[1]-0.5*Step[1];
transform[1] = Step[1];
transform[2] = 0.0;
transform[3] = Start[0]-0.5*Step[0];
transform[4] = 0.0;
transform[5] = Step[0];
Ncont = 5;
/* find median without sorting */
k=0; Max=Min; MedTest = 0.5*(_CHI_SCANS_*_CHI_SCANS_);
for (i = 0; i < _CHI_SCANS_; ++i) {
for (j = 0; j < _CHI_SCANS_; ++j) {
if ((Y[i][j] > Max) && (k <= (int)(MedTest)) ) {
Max = Y[i][j]; ++k;
}
}
}
for (i = 0; i<Ncont; ++i) {
cont_levels[i] = Min + (i+1)*(Max-Min)/(4.0*Ncont);
if (cont_levels[i] <= 9.9999)
sprintf(cont_labels[i], "%6.4f", cont_levels[i]);
}
if (Flags.eps == OFF) {
if(cpgopen("/XSERVE") != 1)
#ifdef _WITH_GTK
{
if (Flags.interactive == ON) {
make_dialog (_(errmsg[2]));
return(8);
} else nf_error(_(errmsg[2]));
}
#else
nf_error(_(errmsg[2]));
#endif
} else {
if(cpgopen(file_out) != 1)
#ifdef _WITH_GTK
{
if (Flags.interactive == ON) {
make_dialog (_(errmsg[1]));
return(8);
} else nf_error(_(errmsg[1]));
}
#else
nf_error(_(errmsg[1]));
#endif
++Flags.eps;
}
#ifdef _WITH_GNUPLOT
gnu_start();
#endif
cpgscf(2); cpgslw(1); cpgsch(1.6);
cpgsvp(0.15, 0.9, 0.15, 0.9);
cpgswin(Start[1], (Start[1] + _CHI_SCANS_ * Step[1]),
Start[0], (Start[0] + _CHI_SCANS_ * Step[0]) );
cpglab( Param[sec_p].Name, Param[fir_p].Name, "");
cpgbox("BCNST", 0,0, "BCNST", 0,0);
cpgcont(a_ptr, _CHI_SCANS_, _CHI_SCANS_, 1, _CHI_SCANS_, 1, _CHI_SCANS_,
cont_levels, (-Ncont), transform);
#ifndef _WITH_GNUPLOT
cpgsch(0.99);
for (i = 0; i < Ncont; ++i)
cpgconl(a_ptr, _CHI_SCANS_, _CHI_SCANS_, 1, _CHI_SCANS_, 1, _CHI_SCANS_,
cont_levels[i], transform,
(char *) &cont_labels[i], (int) 200, (int) 5);
cpgsch(1.6);
#endif
if (Flags.eps != OFF) my_cpgend();
else cpgend();
#endif
return(0);
}
syntax highlighted by Code2HTML, v. 0.9.1