/*      FILE: snowman.c
 *   PURPOSE: 
 *    AUTHOR: Kenneth J. Mighell (mighell@noao.edu)
 *  LANGUAGE: ANSI C
 *      DATE: 2001SEP10
 * COPYRIGHT: (C) 2001 Assoc. of Universities for Research in Astronomy Inc.
 */

#include <mx.h>
#include "inc.h"

#define DEBUG (MX_TRUE)

IRAFTASK (snowman) {
 int i;
 extern int iraf_snowman (void);
 i = iraf_snowman ();
}

/* 
 * ############################################################################
 */
int iraf_snowman(void) {

  char mxfunc[] = "iraf_snowman";
  int mxstatus=0;
  int status;

  inc_ParS ParS;
  inc_ImageS ImageS;
  struct mxip_image_s *xd=NULL;
  struct mxip_image_s *yd=NULL;
  struct mxip_image_s *zd=NULL;
  struct mxip_image_s *ed=NULL;
  struct mxip_image_s *md=NULL;
  struct mxip_image_s *wd=NULL;
  struct mxip_image_s *del[MXNPARL];
  double mnpeak;
  double mxpeak;
  double rond;
  double pixeld;
  long nparl;
  long xdiml;
  long ydiml;
  long nptsl;
  long jl;
  int j;
  char ausec[MXNPARL+1];
  double trued[MXNPARL];
  double guessd[MXNPARL];
  double ad[MXNPARL];
  double ead[MXNPARL];
  double diff_stopd;
  double lambdad;
  double old_chisqd;
  double chisqd;
  double rchsqd;
  double erchsqd;
  int k;
  long nptsmxl;
  double diff;
  long ndfl;
  int mfiti;
  double nsr;

  mxstatus++;
  status = mxp_init_f0();
  if (status) goto mx_error;

  mxstatus++;
  status = inc_ParS_Get_i1 (&ParS);
  if (status) goto mx_error;

  mxstatus++;
  status = inc_ImageS_Open_i2 (&ParS, &ImageS);
  if (status) goto mx_error;

  mxstatus++;
  status = mxip_image_s_open_f7( &xd, 
    ImageS.nxi, ImageS.nyi, ImageS.llx0i, ImageS.lly0i, 1, "X values");
  if (status) goto mx_error;

  mxstatus++;
  status = mxip_image_s_open_f7( &yd, 
    ImageS.nxi, ImageS.nyi, ImageS.llx0i, ImageS.lly0i, 1, "Y values");
  if (status) goto mx_error;

  zd = ImageS.mxIrafImage->data;

  mxstatus++;
  status = mxip_image_s_open_f7( &ed, 
    ImageS.nxi, ImageS.nyi, ImageS.llx0i, ImageS.lly0i, 1, "Error values");
  if (status) goto mx_error;

  mxstatus++;
  status = mxip_image_s_open_f7( &md, 
    ImageS.nxi, ImageS.nyi, ImageS.llx0i, ImageS.lly0i, 1, "Model values");
  if (status) goto mx_error;

  mxstatus++;
  status = mxip_image_s_open_f7( &wd, 
    ImageS.nxi, ImageS.nyi, ImageS.llx0i, ImageS.lly0i, 1, "Working vector");
  if (status) goto mx_error;

  mxstatus++;
  for (k=0; k<MXNPARL; k++) {
    sprintf( MX.tmpmsg, "del[%d]", k);
    status = mxip_image_s_open_f7( &del[k], 
      ImageS.nxi, ImageS.nyi, ImageS.llx0i, ImageS.lly0i, 1, MX.tmpmsg);
    if (status) {
      sprintf( MX.tmpmsg, "Had problems creating del[%d]", k);
      goto mx_error2;
    }
  }
  
  xdiml = ImageS.nxi;
  ydiml = ImageS.nyi;
  nptsl = xdiml*ydiml;
  nptsmxl = nptsl;
  inc_SetXY_v4( xd->vectord, yd->vectord, xdiml, ydiml);
  rond = ParS.ron;
  if (rond<=0) rond = 1e-4;
  for (jl=0; jl<nptsl; jl++) {
    pixeld = zd->vectord[jl];
    if (pixeld>0) {
      ed->vectord[jl] = sqrt(pixeld+(rond*rond)); 
    } else {
      ed->vectord[jl] = rond;
    }
  }
  strcpy(ausec,ParS.parflags);
  for (j=0; j<MXNPARL; j++) {
    if (ausec[j]=='c') guessd[j] = trued[j];
  }
  trued[0] = ParS.tbe;
  trued[1] = ParS.tbx;
  trued[2] = ParS.tby;
  trued[3] = ParS.tse;
  trued[4] = ParS.tsx - 0.5;
  trued[5] = ParS.tsy - 0.5;
  trued[6] = ParS.tss;
  guessd[0] = ParS.gbe;
  guessd[1] = ParS.gbx;
  guessd[2] = ParS.gby;
  guessd[3] = ParS.gse;
  guessd[4] = ParS.gsx-ImageS.llx0i-0.5;
  guessd[5] = ParS.gsy-ImageS.lly0i-0.5;
  guessd[6] = ParS.gss;
  for (j=0; j<MXNPARL; j++) {
    if (ausec[j]=='c') guessd[j] = trued[j];
  }
  for (j=0; j<MXNPARL; j++) {
    ad[j] = guessd[j];
  }
  nparl=7;
  diff_stopd = 0.5;
  lambdad = 0.001;   
  old_chisqd = 1e20;
  for (k=1; k<=10; k++) {
  
    status = LevenbergMarquardt_i18(
      nptsl,
      xd,
      yd,
      zd,
      ed,
      md,
      wd,
      del,
      nparl,
      ad,   
      ead,
      ausec,
      diff_stopd,
      &chisqd,
      &ndfl,
      &mfiti,
      &lambdad,
      ParS.method_is_digital
    );
    if (status) goto mx_error;
    diff = old_chisqd - chisqd;
    if (diff<diff_stopd) goto stop;
    old_chisqd = chisqd;
  }



stop:
  ParS.oniter = k;
  ParS.oebe = ead[0];
  ParS.oebx = ead[1];
  ParS.oeby = ead[2];
  ParS.oese = ead[3];
  ParS.oesx = ead[4];
  ParS.oesy = ead[5];
  ParS.oesf = ead[6]*SIGMA_PER_FWHM;
  ParS.obe = ad[0];
  ParS.obx = ad[1];
  ParS.oby = ad[2];
  ParS.ose = ad[3];
  ParS.osx = ad[4] + ImageS.llx0i + 0.5;
  ParS.osy = ad[5] + ImageS.lly0i + 0.5;
  ParS.osf = ad[6]*SIGMA_PER_FWHM;
  ParS.ochisq = chisqd;
  ParS.ondf = ndfl;
  ParS.orchsq = chisqd/ndfl;
  ParS.oerchsq = sqrt(2.0*ndfl)/ndfl;
  ParS.onpar = nparl;
  ParS.omfit = mfiti;
  ParS.osmag = 44.4444;
  ParS.oesmag = 44.4444;
  if (ParS.ose>0) ParS.osmag = -2.50*log10(ParS.ose/1.0);
  nsr = ParS.oese/ParS.ose;
  if (nsr>0.0) ParS.oesmag = 1.0857*nsr; /* 1.0857 is Pogson's constant */

  /*
  ** =========================================================================
  */

  mxstatus++;
  status = inc_ParS_Put_i1 (&ParS);
  if (status) goto mx_error;

  inc_ParS_Print_v2 (&ParS, stdout);

ok:
  status = 0;
  goto bye;

mx_error:
mx_error1:
  mxp_tmpmsg_init_f0();
mx_error2:
  mxp_errmsg_append_f3 (mxfunc, mxstatus, MX.tmpmsg);
  mxp_tmpmsg_init_f0();
  mxp_errmsg_dump_f0();
  goto bye;

bye:
  inc_ImageS_Close_i1( &ImageS );
  fflush (stdout);
  return(status);

}
#undef DEBUG
/* end-of-file */