/*
Routine to compute mosaic weigthed sky images.
  Sept 2003 Jérôme Chenevez
*/

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<string.h>
#include<jemx.h>
#include<isdc.h>
#include<dal3aux.h>

#define NB_IMAGES 999       /* Maximum number of input images here (may be increased) */
#define XD 600              /* X dim. of input images  (NL: 255) */
#define YD 600              /* Y dim. of input images  (CBJ: New:525, Old:262) */
#define MOSAICSIZE 1000000

/* Declare big arrays as static ! */
  int mark[MOSAICSIZE];
  double pval_o[MOSAICSIZE], var_o[MOSAICSIZE], exp_o[MOSAICSIZE], var_s[MOSAICSIZE], exp_s[MOSAICSIZE], t_exp[MOSAICSIZE], weight[MOSAICSIZE], weigh2[MOSAICSIZE];
  double pvalue[XD*YD], variance[XD*YD];
  double pval_i[XD][YD], var_i[XD][YD], exp_i[XD][YD];
  double vign[XD][YD], vigntmp[XD][YD];
  float value_in[XD*YD], value_out[XD*YD];
  float pvalcor_f[525*525];  /* CBJ's new correction array */

/* ********************************************************** */
#define      sky_xdim     255
#define      sky_ydim     255
#define      skydim       sky_xdim*sky_ydim

  void clean(float sky_r[skydim], float skyim[skydim], int bol);
/* ********************************************************** */

main ()
{
  int i=0, j=0, k=0, loop=0;
  int xdim=XD, ydim=YD, dimx=(int)(sqrt(MOSAICSIZE)), dimy=(int)(sqrt(MOSAICSIZE));
  int NL=0, CBJ=0;
  int numHDU=2 , numAxes=2;
  int p=0, q=0, pix[2]={0, 0};
  int index1=0, index2=0;
  int lin=0, status=ISDC_OK;
  int isub, nsub, psub=0, qsub=0;
  int ext[NB_IMAGES];
  int numImages = 0;  /* Number of images found  */
  int pvalcor[262*262];  /* CBJ's correction array */
  int nex=1;
  long startValues[2] = {1L, 1L};
  long endValues[2] = {(long)(XD), (long)(YD)};
  long numValues=0L;
  double pi, dr;
  float dsub, reduc, dp=0., dq=0.;
  float wei=1., radius=5., diamout=15.;
  double ra_o=0., de_o=0., ra_x=0., de_x=0.;
  double t1[3][3], t2[3][3], t[3][3];
  double t1_o[3][3], t2_o[3][3], t_o[3][3], o_t[3][3];
  double NaN = NAN;
  double ra=0., dec=0.;
  double centerA=-1.0, centerD=0.0;
  double exp_t[NB_IMAGES];
  double crpix1[NB_IMAGES], crpix2[NB_IMAGES];
  double cd1_1[NB_IMAGES], cd2_1[NB_IMAGES],cd1_2[NB_IMAGES], cd2_2[NB_IMAGES];
  double crval1[NB_IMAGES], crval2[NB_IMAGES];
  double crval1_min=0., crval1_max=0., crval2_min=0., crval2_max=0., crval1_o=0., crval2_o=0.;
  double crpix1_o=0., crpix2_o=0.;
  double cdelt=0., cd1_1_o=0., cd1_2_o=0., cd2_1_o=0., cd2_2_o=0.;
  double yz[2], ayz[2], vi[3], ve[3], vi_o[3];
  double maxvign=0.;
  char DOLstring[NB_IMAGES][99], DOLstring_1[NB_IMAGES][99], DOLstring_2[99];  /* DAL_FILENAME_STRING */
  char fitstring[99];
  char parameter[99];
  char mosaic[] = "MOSAIC";
  char vignef[] = "VIGNET";
  char template[] = "JMX2-SKY.-IMA.tpl";
  char skymaps[99];
  char vign_s[99];
  char filexpstring[NB_IMAGES][99];
  char useA='N', plot_T='N', plot_I='N', plot_E='N', plot_V='N', plot_W='N', plot_W2='N', plot_SoE='N', plot_RVE='N';
  void *arrayPtr = NULL;
  void *arrayBuff = NULL;
  FILE *vignet;
  FILE *fp1;

  dal_dataType DALtype = DAL_DOUBLE;
  dal_element *dolPtr = NULL;
  dal_element *skyPtr = NULL;
  dal_element *pvlPtr = NULL;
  dal_element *expPtr = NULL;
  dal_element *varPtr = NULL;
  dal_element *sw_Ptr = NULL;
  dal_element *sw2Ptr = NULL;
  dal_element *mapPtr = NULL;
  dal_element *sqtPtr = NULL;
  dal_element *texPtr = NULL;
  fitsfile *fitsinPtr = NULL;


/* Statements begin */

  pi = acos (-1.);
  dr = pi/180.;
  isub = 2;
  nsub = 2*isub + 1;
  dsub = 1./nsub;
  reduc = dsub*dsub;

/* Read input parameters */
  if( (fp1 = fopen("mosaic_weight.par", "rt")) == NULL ) {
     printf("*** Could not open the necessary parameters file named mosaic_weight.par STOP! ***\n");
     return;
  }
  for (i=0; i<100; i++){
    fgets(parameter, 99, fp1);
    if (feof(fp1)) goto x_end;
    sscanf(parameter, "%d", &lin);
    if (lin == 1) sscanf(parameter, "%d %s", &lin, &skymaps);
    if (lin == 2) sscanf(parameter, "%d %f", &lin, &radius);
    if (lin == 3) sscanf(parameter, "%d %f", &lin, &diamout);
    if (lin == 31) sscanf(parameter, "%d %lf", &lin, &centerA);
    if (lin == 32) sscanf(parameter, "%d %lf", &lin, &centerD);
    if (lin == 101) sscanf(parameter, "%d %c", &lin, &useA);
    if (lin==201) sscanf(parameter, "%d %c", &lin, &plot_T);
    if (lin==202) sscanf(parameter, "%d %c", &lin, &plot_I);
    if (lin==203) sscanf(parameter, "%d %c", &lin, &plot_E);
    if (lin==204) sscanf(parameter, "%d %c", &lin, &plot_V);
    if (lin==205) sscanf(parameter, "%d %c", &lin, &plot_W);
    if (lin==206) sscanf(parameter, "%d %c", &lin, &plot_W2);
    if (lin==207) sscanf(parameter, "%d %c", &lin, &plot_SoE);
    if (lin==208) sscanf(parameter, "%d %c", &lin, &plot_RVE);
  }
x_end:
  fclose (fp1);

/* Look for number and size of input images */
  if( (fp1 = fopen(skymaps, "rt")) == NULL ) {
     printf("Could not open skymaps list file\n");
     return;
  }
  if(1 != fscanf(fp1, "%d", &numImages)){
    fprintf(stderr, "The skymaps list file must start with the number of input fits files to be read! \n");
    return;
  }
    if(1 != fscanf(fp1, "%s", fitstring)){
      fprintf(stderr, "Can't read fits file name in skymaps list file\n");
      return 1;
    }
  strcat(fitstring, "[1]");
  fits_open_file(&fitsinPtr, fitstring, READONLY, &status); /* Open input fits file to get Ptr */
  fits_read_key_lng(fitsinPtr, "NAXIS1", &numValues, NULL, &status); /* Get total number of extensions in the input fits */
  fits_close_file(fitsinPtr, &status); /* Close input fits file */
  rewind(fp1);
  fscanf(fp1, "%d", &numImages);

  if (numImages > NB_IMAGES){
    printf("Please set the parameter NB_IMAGES > %d and recompile the program\n", numImages);
    return;
  }
  else printf("There are %d input images of dimensions %ld * %ld pixels\n", numImages, numValues, numValues);

/* Initialisations */

  ydim = xdim = (int)(numValues);   /* Assuming that the input skymaps have always square dimensions */
  if (xdim == 255) NL = 1;
  else if ((xdim == 262) || (xdim==525)) CBJ = 1;

  endValues[0] = numValues;
  endValues[1] = numValues;

  for (i=0; i<numImages; i++){
    ext[i] = 0;
    crpix1[i] = crpix2[i] = 0;
    exp_t[i] = 0.;
    cd1_1[i] = cd2_1[i] = cd1_2[i] = cd2_2[i] = 0.;
    crval1[i] = crval2[i] = 0.;
  }
  for (i=0; i<MOSAICSIZE; i++){
    pval_o[i] = 0.;
    var_o[i] = 0.;
    exp_o[i] = 0.;
    var_s[i] = 0.;
    exp_s[i] = 0.;
    t_exp[i] = 0.;
    weight[i] = 0.;
    weigh2[i] = 0.;
  }
  for (i=0; i<xdim; i++)
    for (j=0; j<ydim; j++){
      pval_i[i][j] = NaN;
      var_i[i][j] = 0.;
      exp_i[i][j] = 0.;
      vign[i][j] = 0.;
    }
  for (i=0; i<xdim*ydim; i++){
    pvalue[i] = NaN;
    variance[i] = 0.;
  }
  for (i=0; i<3; i++){
    vi_o[i] = 0.;
    ve[i] = 0.;
    for (j=0; j<3; j++){
      t1[i][j] = t2[i][j] = t[i][j] =0.;
      t1_o[i][j] = t2_o[i][j] = t_o[i][j] = 0.;
      o_t[i][j] = 0.;
    }
  }

/* Initialisations done */

/* Read vignetting-function array (in mm2) from txt file (NL or CBJ) */
  if (NL == 1) strcpy(vign_s, "illum_k-2.txt");
  else if (CBJ == 1){        /* Read New vignetting factors from CBJ */
    if (xdim == 262){
      if( (vignet = fopen("cbj_vig.dat","rt")) == NULL ) {
        printf("Could not open the file named: cbj_vig.dat\n");
        return(0);
      }
      j = 0;
      while (j<xdim*ydim){
        for (i=0; i<43; i++){
          fgets(vign_s, 99, vignet);
          sscanf(vign_s, "%d %d %d %d %d %d", &pvalcor[j++],&pvalcor[j++],&pvalcor[j++],&pvalcor[j++],&pvalcor[j++],&pvalcor[j++]);
        }
        fgets(vign_s, 99, vignet);
        sscanf(vign_s, "%d %d %d %d", &pvalcor[j++],&pvalcor[j++],&pvalcor[j++],&pvalcor[j++]);
      }
      fclose(vignet);
      strcpy(vign_s, "vignet_262_2+.txt");
    }
    else if (xdim == 525){
      if( (vignet = fopen("bckim_1_5.txt","rt")) == NULL ) {
        printf("Could not open the file named: bckim_1_5.txt\n");
        return(0);
      }
      j = 0;
      while (j<xdim*ydim){
        for (i=0; i<131; i++){
          fgets(vign_s, 99, vignet);
          sscanf(vign_s, "%f %f %f %f", &pvalcor_f[j++],&pvalcor_f[j++],&pvalcor_f[j++],&pvalcor_f[j++]);
        }
        fgets(vign_s, 99, vignet);
        sscanf(vign_s, "%f", &pvalcor_f[j++]);
      }
      fclose(vignet);
      strcpy(vign_s, "vignet_525_2+.txt");
    }
  }
  printf("Use of %s for vignetting array\n", vign_s);
  if( (vignet = fopen(vign_s,"rt")) == NULL ) {
    printf("Could not open %s\n", vign_s);
    return(0);
  }
  for (i=0; i<xdim; i++)
    for (j=0; j<ydim; j++)
      if(1 !=  fscanf(vignet,"%lf", &vign[i][j])){
        fprintf(stderr,"Can't read vignetting at [%d,%d]\n", i, j);
        return 1;
      }
  fclose(vignet);
  /* Smoothing */
  for (i=1; i<xdim-1; i++)
    for (j=1; j<ydim-1; j++){
      vigntmp[i][j] = (vign[i-1][j-1]+vign[i][j-1]+vign[i+1][j-1]+\
      vign[i+1][j]  +vign[i][j]  +vign[i-1][j]+\
      vign[i-1][j+1]+vign[i][j+1]+vign[i+1][j+1])/9.;
      if (vigntmp[i][j] > maxvign) maxvign = vigntmp[i][j];
    }
  for (i=0; i<xdim; i++)                  /* Normalisation */
    for (j=0; j<ydim; j++) vign[i][j] = vigntmp[i][j]/maxvign;
/* Convert vignetting array to FITS format
  status = DALobjectCreate(vignef, DAL_DISK, template, &skyPtr, status);
  if (status != ISDC_OK){
    printf("Cannot create vignet object, status = %d\n",status);
    return status;
  }
  status = DALarrayModStruct(skyPtr, DALtype, numAxes, endValues, status);
  if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
  arrayBuff = vign;
  status = DALarrayPutSection(skyPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
  if (status != ISDC_OK){
    printf("Cannot put section vign object, status = %d\n",status);
    return status;
  }
  status = DALobjectClose(skyPtr, DAL_SAVE, status);
  if (status != ISDC_OK){
    printf("Cannot close object vignef! status = %d\n",status);
    return status;
  }
*/

/* Read HEADER attributes of each input images */
  for (loop=0; loop<numImages; loop++){
  /* Loop over input skymaps */
    if(1 != fscanf(fp1, "%s", DOLstring[loop])){
      fprintf(stderr, "Can't read fits file name in skymaps list file\n");
      return 1;
    }
    if (strcmp(DOLstring[loop], DOLstring[loop-1]) == 0){
      printf("*** Stop: DOL %s is repeated! ***\n", DOLstring[loop]);
      return 1;
    }
    strcpy(DOLstring_1[loop], DOLstring[loop]);

    if (CBJ == 1) strcat( DOLstring_1[loop], "[1]");
    else if (NL == 1) {
      fits_open_file(&fitsinPtr, DOLstring[loop], READONLY, &status); /* Open input fits file to get Ptr */
      fits_get_num_hdus(fitsinPtr, &numHDU, &status); /* Get total number of extensions in the input fits */
      fits_close_file(fitsinPtr, &status); /* Close input fits file */
      if (useA == 'N') strcat( DOLstring_1[loop], "[1]");
      if ((useA == 'Y') && (numHDU <= 2)) strcat( DOLstring_1[loop], "[1]");
      else if (useA =='Y'){
        strcat( DOLstring_1[loop], "[2]");
        ext[loop] = 2;
      }
    }

    printf("The %d. DOL is : %s\n", loop+1, DOLstring_1[loop]);
    status = DALobjectOpen(DOLstring_1[loop], &dolPtr, status);
    if (status != ISDC_OK){
      printf("Cannot open DOL %s for image number %d\n", DOLstring_1[loop], loop+1);
      return status;
    }
    status = DALattributeGetReal(dolPtr, "ONTIME", &exp_t[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CD1_1 for image number %d\n", loop+1);
      return( status );
    }
    if (exp_t[loop] <= 0.) printf("*** Warning: Image %d, ONTIME = %f ***\n", loop+1, exp_t[loop]);
    status = DALattributeGetReal(dolPtr, "CD1_1", &cd1_1[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CD1_1 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CD2_1", &cd2_1[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CD2_1 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CD1_2", &cd1_2[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CD1_2 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CD2_2", &cd2_2[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CD2_2 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CRVAL1", &crval1[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CRVAL1 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CRVAL2", &crval2[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CRVAL2 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CRPIX1", &crpix1[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CRPIX1 for image number %d\n", loop+1);
      return( status );
    }
    status = DALattributeGetReal(dolPtr, "CRPIX2", &crpix2[loop], NULL, NULL, status);
    if(status != ISDC_OK){
      printf("Cannot read attribute CRPIX2 for image number %d\n", loop+1);
      return( status );
    }
    status = DALobjectClose(dolPtr, DAL_SAVE, status);
    if (status != ISDC_OK){
      printf("Cannot close DOL for image number %d\n", loop+1);
      return status;
    }
  }

  fclose (fp1);

/* Define the output mosaic skymap */
  crval1_min = 360.;
  crval1_max = 0.;
  crval2_min = 90.;
  crval2_max = -90.;
  for(loop=0; loop<numImages; loop++){
    if (crval1[loop] < crval1_min) crval1_min = crval1[loop];
    if (crval2[loop] < crval2_min) crval2_min = crval2[loop];
    if (crval1[loop] > crval1_max) crval1_max = crval1[loop];
    if (crval2[loop] > crval2_max) crval2_max = crval2[loop];
  }
  if ((crval1_max-crval1_min) >= 180.) crval1_max -= 360.;
  if (diamout > 0.)
  if ( sqrt((crval1_max-crval1_min)*(crval1_max-crval1_min) + (crval2_max-crval2_min)*(crval2_max-crval2_min)) >= diamout) 
    printf("*** Warning: total separation angle larger than %f degrees, some images may not be included ***\n", diamout);
  if (centerA < 0.){
    crval1_o = 0.5 * (crval1_min + crval1_max);
    if (crval1_o < 0.) crval1_o+= 360.;
    crval2_o = 0.5 * (crval2_min + crval2_max);
  }
  else {
  crval1_o = centerA;
  crval2_o = centerD;
  }
  printf("Center of output mosaic skyimage: RA=%f, DEC=%f\n", crval1_o, crval2_o);

/* Define transformation from sky coord. to FITS */
  cdelt = 0.03 ; /* degrees per pixel*/
  if (diamout > 0.){
    dimx = (int)(diamout / cdelt);
    dimy = (int)(diamout / cdelt);
  }
  crpix1_o = (float)(dimx)/2.;
  crpix2_o = (float)(dimy)/2.;
  cd1_1_o = - cdelt;
  cd1_2_o = 0.0;
  cd2_1_o = 0.0;
  cd2_2_o = cdelt;
  ra_o = crval1_o;
  de_o = crval2_o;
  t1_o[0][0] = cos(ra_o*dr);
  t1_o[1][0] = sin(ra_o*dr);
  t1_o[2][0] = 0.;
  t1_o[0][1] = -sin(ra_o*dr);
  t1_o[1][1] = cos(ra_o*dr);
  t1_o[2][1] = 0.;
  t1_o[0][2] = 0.;
  t1_o[1][2] = 0.;
  t1_o[2][2] = 1.;
  t2_o[0][0] = cos(de_o*dr);
  t2_o[1][0] = 0.;
  t2_o[2][0] = sin(de_o*dr);
  t2_o[0][1] = 0.;
  t2_o[1][1] = 1.;
  t2_o[2][1] = 0.;
  t2_o[0][2] = -sin(de_o*dr);
  t2_o[1][2] = 0.;
  t2_o[2][2] = cos(de_o*dr);
  for (i=0; i<3; i++)
    for (j=0; j<3; j++){
      t_o[i][j] = 0.;
      for (k=0; k<3; k++) t_o[i][j] += t1_o[i][k]*t2_o[k][j];
    }
  for (i=0; i<3; i++)
    for (j=0; j<3; j++) o_t[i][j] = t_o[j][i];

/* COMPOSE MOSAIC BY LOOPING OVER SKY IMAGES TO ADD */
  for (loop=0; loop<numImages; loop++){
    if ( diamout > 0.)
    if ((sqrt((crval1[loop]-crval1_o)*(crval1[loop]-crval1_o)+(crval2[loop]-crval2_o)*(crval2[loop]-crval2_o))-radius) >= diamout/2.){
      printf(" Excluding %s\n", DOLstring[loop]);
      continue;
    }
    for (i=0; i<MOSAICSIZE; i++) mark[i] = 0;
  /* Define transformation from FITS to sky coord. */
    ra_x = crval1[loop];
    de_x = crval2[loop];
    t1[0][0] = cos(ra_x*dr);
    t1[1][0] = sin(ra_x*dr);
    t1[2][0] = 0.;
    t1[0][1] = -sin(ra_x*dr);
    t1[1][1] = cos(ra_x*dr);
    t1[2][1] = 0.;
    t1[0][2] = 0.;
    t1[1][2] = 0.;
    t1[2][2] = 1.;
    t2[0][0] = cos(de_x*dr);
    t2[1][0] = 0.;
    t2[2][0] = sin(de_x*dr);
    t2[0][1] = 0.;
    t2[1][1] = 1.;
    t2[2][1] = 0.;
    t2[0][2] = -sin(de_x*dr);
    t2[1][2] = 0.;
    t2[2][2] = cos(de_x*dr);
    for (i=0; i<3; i++)
      for (j=0; j<3; j++){
        t[i][j] = 0.;
        for (k=0; k<3; k++){
          t[i][j] += t1[i][k]*t2[k][j];
        }
      }

  /* Read skymaps inputs */
    status = DALobjectOpen(DOLstring_1[loop], &dolPtr, status);
    if (status != ISDC_OK){
      printf("Cannot open DOL %s for image number %d\n", DOLstring_1[loop], loop+1);
      return status;
    }
    arrayPtr = pvalue;
    status = DALarrayGetSection(dolPtr, numAxes, startValues, endValues, &DALtype, &numValues, arrayPtr, status);
    if (status != ISDC_OK){
      printf("Cannot get section pvalues for image number %d\n", loop+1);
      return status;
    }
    status = DALobjectClose(dolPtr, DAL_SAVE, status);
    if (status != ISDC_OK){
      printf("Cannot close DOL for image number %d\n", loop+1);
      return status;
    }

    if ((NL == 1) && (ext[loop] == 2)){  /* Use of NL filtering function for large scale artefacts in images without strong sources */
      for (p=0; p<(int)(numValues); p++) value_in[p] = (float)(pvalue[p]);
      arrayPtr = value_out;
      clean(value_in, arrayPtr,1);
      for (p=0; p<(int)(numValues); p++) pvalue[p] = (double)(value_out[p]);
    }

    strcpy(DOLstring_2, DOLstring[loop]);
    if (NL == 1) strcat(DOLstring_2, "[0]");              /* Extension is [0] for NL's variance maps */
    else if (CBJ == 1) strcat(DOLstring_2, "[2]");
    printf("Variance map is in : %s\n", DOLstring_2);
    status = DALobjectOpen(DOLstring_2, &dolPtr, status);
    if (status != ISDC_OK){
      printf("Cannot open DOL %s for image number %d\n", DOLstring_2, loop+1);
      return status;
    }
    arrayPtr = variance;
    status = DALarrayGetSection(dolPtr, numAxes, startValues, endValues, &DALtype, &numValues, arrayPtr, status);
    if (status != ISDC_OK){
      printf("Cannot get section variance for image number %d\n", loop);
      return status;
    }
    status = DALobjectClose(dolPtr, DAL_SAVE, status);
    if (status != ISDC_OK){
      printf("Cannot close DOL for image number %d\n", loop+1);
      return status;
    }
    if (NL == 1){  /* Use of NL filtering function for large scale artefacts in variance maps */
      for (p=0; p<(int)(numValues); p++) value_in[p] = (float)(variance[p]);
      arrayPtr = value_out;
      clean(value_in, arrayPtr,0);
      for (p=0; p<(int)(numValues); p++) variance[p] = (double)(value_out[p]);
    }
     
    for (p=0; p<(int)(numValues); p++){
      pvalue[p] = pvalue[p] * reduc;  /* Divide the pixel values by the nb of pixel subdivisions */
      variance[p] = variance[p] * reduc; /* Divide the variances by the nb of pixel subdivisions */
      if ((CBJ == 1) && (xdim == 262)) pvalue[p] = pvalue[p] * pvalcor[p];
      else if ((CBJ == 1) && (xdim == 525)) pvalue[p] = pvalue[p] * pvalcor_f[p];
    }

    /* Loop over individual pixels */
    for (p=0; p<xdim; p++){
      for (q=0; q<ydim; q++){
        index1 = q*xdim+p;
        pval_i[p][q] = isnan(pvalue[index1]) ? 0.0 : pvalue[index1];
        var_i[p][q] = variance[index1] +100.;  /* Add a threshold to the variance to avoid too high weights */
        exp_i[p][q] = exp_t[loop]*reduc * vign[p][q]; /* Compute effective exposure with vignetting-function in cm2 */
        yz[0] = cd1_1[loop]*(p+1-crpix1[loop]) + cd1_2[loop]*(q+1-crpix2[loop]);
        yz[1] = cd2_1[loop]*(p+1-crpix1[loop]) + cd2_2[loop]*(q+1-crpix2[loop]);
        if (sqrt(yz[0]*yz[0]+yz[1]*yz[1]) > radius) continue;
        /* Loop over pixels subdivisions  */
        for (psub=-isub; psub<=isub; psub++){
        for (qsub=-isub; qsub<=isub; qsub++){
          dp = psub*dsub;
          dq = qsub*dsub;
          yz[0] = cd1_1[loop]*(p+dp+1-crpix1[loop]) + cd1_2[loop]*(q+dq+1-crpix2[loop]);
          yz[1] = cd2_1[loop]*(p+dp+1-crpix1[loop]) + cd2_2[loop]*(q+dq+1-crpix2[loop]);
          for(k=0; k<2; k++) ayz[k] = atan(yz[k]*dr);
          for(k=1; k<3; k++) vi[k] = sin(ayz[k-1]);
          vi[0] = sqrt(1. -vi[1]*vi[1]-vi[2]*vi[2]);
          for (i=0; i<3; i++){
            ve[i] = 0.;
            for (j=0; j<3; j++) ve[i] += t[i][j]*vi[j];
          }
/*
          dec = asin(ve[2])/dr;
          ra = atan2(ve[1],ve[0])/dr + 360.;
          if (ra > 360.) ra -= 360.;
*/
          /* Where does this position belong in mosaic image? */
          for (i=0; i<3; i++){
            vi_o[i] = 0.;
            for (j=0; j<3; j++) vi_o[i] += o_t[i][j]*ve[j];
          }
          for (k=0; k<2; k++){
            ayz[k] = asin(vi_o[k+1]);
            yz[k] = tan(ayz[k]) / dr;
          }
          pix[0] = (int)((-1./cdelt)*yz[0]+ crpix1_o -1.);
          pix[1] = (int)((1./cdelt)*yz[1] + crpix2_o -1.);
          index2 = pix[1]*dimx+pix[0];
          if ((pix[0]>=0) && (pix[0]<dimx) && (pix[1]>=0) && (pix[1]<dimy)){
            if (mark[index2]==0){
              t_exp[index2] += exp_t[loop]; /* Compute raw exposure times */
              mark[index2] ++;
            }
/*      COMPUTE WEIGTHS AND UPDATE MOSAIC SKYMAP 
            if (exp_o[index2] > 0.)
             wei = (exp_i[p][q]/var_i[p][q]) / (exp_o[index2]/var_o[index2]);
            else wei = 1.0;
            if (wei <= 1.0){
              pval_o[index2] += wei * pval_i[p][q];
              exp_o[index2] += wei * exp_i[p][q];
              var_o[index2] += wei * wei * var_i[p][q];
            }
            else{
              wei = 1./wei; 
              pval_o[index2] = pval_i[p][q] + wei * pval_o[index2];
              if ((exp_i[p][q]+wei*exp_o[index2]) > exp_o[index2]) exp_o[index2] = exp_i[p][q] + wei * exp_o[index2];
              var_o[index2] = var_i[p][q] + wei * wei * var_o[index2];
            }
            weight[index2] += wei;
            weigh2[index2] += wei*wei;
*/        wei = (exp_i[p][q]/var_i[p][q]);
          pval_o[index2] += wei * pval_i[p][q];
          exp_o[index2] += wei * exp_i[p][q];
          var_o[index2] += wei * wei * var_i[p][q];
          weight[index2] += wei;
          weigh2[index2] += wei*wei;
          exp_s[index2] += exp_i[p][q];
          var_s[index2] += var_i[p][q];
          }
        }
        }
      }
    }
  } /* End loop over total number of sky images to add */
  for (index2=0; index2<MOSAICSIZE; index2++){
    pval_o[index2] = pval_o[index2]*var_s[index2]/exp_s[index2];
    exp_o[index2] = exp_o[index2]*var_s[index2]/exp_s[index2];
    var_o[index2] = var_o[index2]*var_s[index2]/exp_s[index2];
    weight[index2] = weight[index2]*var_s[index2]/exp_s[index2];
    weigh2[index2] = weigh2[index2]*var_s[index2]/exp_s[index2];
    if (!pval_o[index2]) pval_o[index2] = NaN;
    if (!exp_o[index2]) exp_o[index2] = NaN;
    if (!var_o[index2]) var_o[index2] = NaN;
    if (!t_exp[index2]) t_exp[index2] = NaN;
    if (!weight[index2]) weight[index2] = NaN;
    if (!weigh2[index2]) weigh2[index2] = NaN;
  }

/* Generate output Mosaic fits file */
  endValues[0] = dimx;
  endValues[1] = dimy;
  if (plot_T == 'Y'){
    sprintf(fitstring, "OBS_TIME_%d", numImages);
    status = DALobjectCreate(fitstring, DAL_DISK, template, &texPtr, status);
    if (status != ISDC_OK){
      printf("Cannot create OBS_TIME object, status = %d\n", status);
      return status;
    }
    printf("Plot exposure Time map in %s.fits\n", fitstring);
    status = DALarrayModStruct(texPtr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar(texPtr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar(texPtr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal(texPtr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal(texPtr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal(texPtr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal(texPtr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal(texPtr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal(texPtr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal(texPtr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal(texPtr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar(texPtr, "IMATYPE", "EXP_TIME", "seconds", "The image shows observation times", status);
    arrayBuff = t_exp;
    status = DALarrayPutSection(texPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section t_exp object\n");
      return status;
    }
    status = DALobjectClose(texPtr, DAL_SAVE, status);
    if (status != ISDC_OK){
      printf("Cannot close object!\n");
      return status;
    }
  }

  sprintf(fitstring, "MOSAIC_%d", numImages);
  status = DALobjectCreate(fitstring, DAL_DISK, "JMX2-SKY.-IMA-IDX.tpl", &skyPtr, status);
  if (status != ISDC_OK){
    printf("Cannot create mosaic object, status = %d\n", status);
    return status;
  }

  if (plot_I == 'Y'){
    nex++;
    printf("Plot raw Intensity map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &pvlPtr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for pval_o object\n");
      return status;
    }
    status = DALarrayModStruct( pvlPtr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( pvlPtr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( pvlPtr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( pvlPtr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( pvlPtr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( pvlPtr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( pvlPtr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( pvlPtr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( pvlPtr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( pvlPtr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( pvlPtr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( pvlPtr, "IMATYPE", "INTENSITY", " ", "The image shows intensity", status);
    arrayBuff = pval_o;
    status = DALarrayPutSection(pvlPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section pval_o object\n");
      return status;
    }
  }
  
  if (plot_E == 'Y'){
    nex++;
    printf("Plot effective Exposure map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &expPtr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for exp_o object\n");
      return status;
    }
    status = DALarrayModStruct( expPtr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( expPtr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( expPtr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( expPtr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( expPtr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( expPtr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( expPtr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( expPtr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( expPtr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( expPtr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( expPtr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( expPtr, "IMATYPE", "EXPOSURE", " ", "The image shows effective exposure", status);
    arrayBuff = exp_o;
    status = DALarrayPutSection(expPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section exp_o object\n");
      return status;
    }
  }

  if (plot_V == 'Y'){
    nex++;
    printf("Plot Variance map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &varPtr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for var_o object\n");
      return status;
    }
    status = DALarrayModStruct( varPtr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( varPtr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( varPtr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( varPtr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( varPtr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( varPtr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( varPtr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( varPtr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( varPtr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( varPtr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( varPtr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( varPtr, "IMATYPE", "VARIANCE", " ", "The image shows variance", status);
    arrayBuff = var_o;
    status = DALarrayPutSection(varPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section var_o object\n");
      return status;
    }
  }

  if (plot_W == 'Y'){
    nex++;
    printf("Plot sum of Weights map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &sw_Ptr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for weight object\n");
      return status;
    }
    status = DALarrayModStruct( sw_Ptr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( sw_Ptr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( sw_Ptr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( sw_Ptr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( sw_Ptr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( sw_Ptr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw_Ptr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw_Ptr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw_Ptr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw_Ptr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( sw_Ptr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( sw_Ptr, "IMATYPE", "WEIGHTS", " ", "The image shows sum of weights", status);
    arrayBuff = weight;
    status = DALarrayPutSection(sw_Ptr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section weight object\n");
      return status;
    }
  }

  if (plot_W2 == 'Y'){
    nex++;
    printf("Plot sum of weights^2 map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &sw2Ptr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for weigh2 object\n");
      return status;
    }
    status = DALarrayModStruct( sw2Ptr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( sw2Ptr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( sw2Ptr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( sw2Ptr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( sw2Ptr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( sw2Ptr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw2Ptr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw2Ptr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw2Ptr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( sw2Ptr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( sw2Ptr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( sw2Ptr, "IMATYPE", "WEIGHT2", " ", "The image shows sum of weights^2", status);
    arrayBuff = weigh2;
    status = DALarrayPutSection(sw2Ptr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section weigh2 object\n");
      return status;
    }
  }

  for (index2=0; index2<MOSAICSIZE; index2++)
    if (exp_o[index2] != 0.){
      pval_o[index2] = pval_o[index2]/exp_o[index2];              /* Actual quantities to plot */
      var_o[index2] = sqrt(var_o[index2])/exp_o[index2];          /* pv/exp and  sqrt(var)/exp */
    }

  if (plot_SoE == 'Y'){
    nex++;
    printf("Plot Signal/Exposure map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &mapPtr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for var_o object\n");
      return status;
    }
    status = DALarrayModStruct( mapPtr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( mapPtr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( mapPtr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( mapPtr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( mapPtr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( mapPtr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( mapPtr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( mapPtr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( mapPtr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( mapPtr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( mapPtr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( mapPtr, "IMATYPE", "SIGNAL", " ", "The image shows signal/exposure", status);
    arrayBuff = pval_o;
    status = DALarrayPutSection(mapPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section new pval_o object\n");
      return status;
    }
  }

  if (plot_RVE == 'Y'){
    nex++;
    printf("Plot SQRT(Variance)/Exposure map in %s.fits[%d]\n", fitstring, nex);
    status = DAL3GENindexCreateMember(skyPtr, NULL, NULL, &sqtPtr, status);
    if (status != ISDC_OK){
      printf("DAL3GEN: Cannot create member for var_o object\n");
      return status;
    }
    status = DALarrayModStruct( sqtPtr, DALtype, numAxes, endValues, status);
    if( status != ISDC_OK )  printf("DALarrayModStruct status : %d\n", status);
    status += DALattributePutChar( sqtPtr, "CTYPE1", "RA---TAN","1962","Projection name", status);
    status += DALattributePutChar( sqtPtr, "CTYPE2", "DEC--TAN","1962 ","Projection name", status);
    status += DALattributePutReal( sqtPtr, "CRVAL1", crval1_o,"deg","Sky coord. ref. 1. axis", status);
    status += DALattributePutReal( sqtPtr, "CRVAL2", crval2_o,"deg","Sky coord. ref. 2. axis", status);
    status += DALattributePutReal( sqtPtr, "CD1_1", cd1_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( sqtPtr, "CD2_2", cd2_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( sqtPtr, "CD1_2", cd1_2_o,"deg", "Transformation", status);
    status += DALattributePutReal( sqtPtr, "CD2_1", cd2_1_o,"deg", "Transformation", status);
    status += DALattributePutReal( sqtPtr, "CRPIX1", crpix1_o,"pixel", "Sky pixel ref. 1. axis", status);
    status += DALattributePutReal( sqtPtr, "CRPIX2", crpix2_o,"pixel", "Sky pixel ref. 2. axis", status);
    status += DALattributePutChar( sqtPtr, "IMATYPE", "SQTVAR", " ", "The image shows sqt(var)/exposure", status);
    arrayBuff = var_o;
    status = DALarrayPutSection(sqtPtr, numAxes, startValues, endValues, DALtype, arrayBuff, status);
    if (status != ISDC_OK){
      printf("Cannot put section new var_o object\n");
      return status;
    }
  }

  status = DALobjectClose(skyPtr, DAL_SAVE, status);
  if (status != ISDC_OK){
    printf("Cannot close object!\n");
    return status;
  }

  return 0;
}

/*--------------------------------------------------------------------------------------*/

void clean(float sky_r[skydim], float skyim[skydim], int bol)
/* this subroutine removes large scale structure from a skymap
   The input skymap is provided in the in the array sky_r
   The cleaned skymap is returned in the array skyim
   Niels Lund, September 2003
*/
{
   int i, j, n, k, max_j[3] = {-1, -1, -1}, min_j[2] = {100000, 100000};
   int indx[10] = {5, sky_ydim-4, sky_ydim+4, 2*sky_ydim-3, 2*sky_ydim+3,
                  3*sky_ydim-2, 3*sky_ydim+2, 4*sky_ydim-1, 4*sky_ydim+1, 5*sky_ydim};
   float b, v, max[2] = {-100000.0, -100000.0}, min[2] = {100000.0, 100000.0};
   
   for (i=5*sky_ydim; i<skydim-5*sky_ydim; i++) {
      for (j=0; j<10; j++) {
         v = sky_r[i+indx[j]];
         if (v > max[0]) {
            max[1] = max[0];
            max_j[1] = max_j[0];
            max[0] = v;
            max_j[0] = i+indx[j];
         } else {
            if (v > max[1]) {
               max[1] = v;
               max_j[1] = i+indx[j];
            }
         }
         if (v < min[0]) {
            min[1] = min[0];
            min_j[1] = min_j[0];
            min[0] = v;
            min_j[0] = i+indx[j];
         } else {
            if (v < min[1]) {
               min[1] = v;
               min_j[1] = i+indx[j];
            }
         }

         v = sky_r[i-indx[j]];
         if (v > max[0]) {
            max[1] = max[0];
            max_j[1] = max_j[0];
            max[0] = v;
            max_j[0] = i-indx[j];
         } else {
            if (v > max[1]) {
               max[1] = v;
               max_j[1] = i-indx[j];
            }
         }
         if (v < min[0]) {
            min[1] = min[0];
            min_j[1] = min_j[0];
            min[0] = v;
            min_j[0] = i+indx[j];
         } else {
            if (v < min[1]) {
               min[1] = v;
               min_j[1] = i+indx[j];
            }
         }
      }
   
      b = 0.0;
      n = 0;
      for (j=0; j<10; j++) {
         k = i+indx[j];
         if (k == max_j[0]) goto try2;
         if (k == max_j[1]) goto try2;
         if (k == min_j[0]) goto try2;
         if (k == min_j[1]) goto try2;
         b += sky_r[k];
         n++;
try2:    k = i-indx[j];
         if (k == max_j[0]) continue;
         if (k == max_j[1]) continue;
         if (k == min_j[0]) continue;
         if (k == min_j[1]) continue;
         b += sky_r[k];
         n++;
      }
      
      b /= (float)(n);
      
      if (bol == 1) skyim[i] = sky_r[i] - b;
      else if (bol == 0) skyim[i] = b;
      else printf("ERROR: bol = %d\n", bol);
   }
   
   return;
}

/* ********************************************************** */