/* 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);

}