#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>
#include <sys/resource.h>
#include <fftw3.h>
#include "jsoc_main.h"
#include "fresize.h"

#define kNOTSPECIFIED "not specified"
#define ARRLENGTH(ARR)  (sizeof(ARR)/sizeof(ARR[0]))
#define PI              (M_PI)
#define RADSINDEG 	(PI/180)
#define RTRUE		(6.96000000e8)  // meters
#define AU		(149597870691)  // meters/au
#define TAU_A           (499.004783806) // light travel time in seconds, = 1 AU/c
#define C               (299792458)     // speed of light, m/s

char *module_name = "mkylms";
char *cvsinfo_mkylms = "cvsinfo: $Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/mkylms.c,v 1.5 2013/05/03 20:44:38 tplarson Exp $";

ModuleArgs_t module_args[] = 
{
   {ARG_STRING, "out",          NULL,           "output data series"},
   {ARG_STRING, "histlink",     "HISTORY",      "name of link to ancillary dataseries for processing metadata"},
   {ARG_INT,    "PERM",         "1",            "set to 0 for transient records, nonzero for permanent records"},
   {ARG_INT,    "VERB",         "1",            "option for level of verbosity: 0=only error and warning messages; >0=print messages outside of loop; >1=print messages inside loop; >2=debugging output", NULL},
   {ARG_STRING, "MODELIST", kNOTSPECIFIED,      "file specifying modes to generate"},
   {ARG_INT,    "XPIXELS",     "1024",          "number of pixels in x direction"},
   {ARG_INT,    "YPIXELS",     "1024",          "number of pixels in y direction"},
   {ARG_FLOAT,  "SOLARP",      "0.0",           "p-angle in degrees"},
   {ARG_FLOAT,  "OBSB0",       "0.0",           "b-angle in degrees"},
   {ARG_FLOAT,  "OBSDIST",     "1.0",           "distance from the sun in au"},
   {ARG_FLOAT,  "XOFFSET",     "0.0",           "offset in pixels from center in x direction"},
   {ARG_FLOAT,  "YOFFSET",     "0.0",           "offset in pixels from center in y direction"},
   {ARG_INT,    "IINTEN",      "0",             "flag for making intensity images, make velocity images otherwise"},
   {ARG_FLOAT,  "PHASE",       "0.0",           "phase in radians"},
   {ARG_FLOAT,  "FREQUENCY",   "0.0",           "frequency in Hz"},
   {ARG_INT,    "ILIMBDARK",   "0",             "flag to correct for limb darkening in intensity images"},
   {ARG_DOUBLES,"coefs",       "0.0",           "limb darkening coeficients, 6 needed"},
   {ARG_INT,    "IHEIGHT",     "0",             "flag to correct for height of formation"},

   {ARG_FLOAT, "AIRY_CDOWN",  "0.0", "ratio of the pixel nyquist frequency to the cutoff frequency of the Airy function"},
   {ARG_INT,   "AIRY_IAP",   "0",    "option for type of apodization"},
   {ARG_INT,   "AIRY_NAP",   "1",    "distance between sampled points in pixels"},
   {ARG_INT,   "AIRY_WID",   "-1",   "half width of kernel"},
   {ARG_INT,   "AIRY_NSUB",  "1",    "distance between sampled points in pixels"},
   {ARG_INT,   "AIRY_XOFF",  "0",    "offset in pixels to add to x index of input image"},
   {ARG_INT,   "AIRY_YOFF",  "0",    "offset in pixels to add to y index of input image"},

   {ARG_INT,   "NBIN",        "-1",   "factor for binning"},
   {ARG_INT,   "BIN_XOFF",    "0",    "offset in pixels in x direction to apply before binning"},
   {ARG_INT,   "BIN_YOFF",    "0",    "offset in pixels in y direction to apply before binning"},

   {ARG_FLOAT, "GAUSS_SIG",   "-1.0", "width of gaussian"},
   {ARG_INT,   "GAUSS_WID",   "-1",   "half width of kernel"},
   {ARG_INT,   "GAUSS_NSUB",  "1",    "distance between sampled points in pixels"},
   {ARG_INT,   "GAUSS_XOFF",  "0",    "offset in pixels to add to x index of input image"},
   {ARG_INT,   "GAUSS_YOFF",  "0",    "offset in pixels to add to y index of input image"},

/*
   {ARG_INT,    "IHORIZ",      "0",             "flag for outputting horizontal component only, otherwise output radial component only"},
   {ARG_INT,    "IBOX",        "0",             "flag for convolving psf with box"},
   {ARG_INT,    "IGAUSS",      "0",             "flag for convolving psf with gaussian"},
   {ARG_FLOAT,  "WPSFX",       "1.0",           "width of gaussian psf in x direction"},
   {ARG_FLOAT,  "WPSFY",       "1.0",           "width of gaussian psf in y direction"},
   {ARG_INT,    "ILORENTZ",    "0",             "flag for convolving psf with lorentzian"},
   {ARG_FLOAT,  "WLORENTZ",    "1.0",           "width of lorenztian psf"},
   {ARG_INT,    "ILININT",     "0",             "flag for convolving psf with linear interpolation"},
   {ARG_INT,    "ICUBINT",     "0",             "flag for convolving psf with cubic interpolation"},
   {ARG_INT,    "IBIN",        "0",             "flag for binning"},
   {ARG_INT,    "NBIN",        "1",             "factor for binning"},
   {ARG_INT,    "ISKIP",       "0",             "flag for subsampling"},
   {ARG_INT,    "NOFF",        "2",             "offset for subsampling"},
*/
   {ARG_TIME,   "TSTART", "1996.05.26_16:00:00_TAI", "first output time"},
   {ARG_INT,    "DTOFF",       "0",             "number of timesteps to offset first image"},
   {ARG_INT,    "NDT",         "0",             "number of time points to generate"},

   {ARG_END}
};

#include "saveparm.c"
#include "timing.c"
#include "set_history.c"
#include "setplm2.c"

static void Ccker(double *u, double s);
static double Ccintd(double *f, int nx, double x);

int DoIt(void)
{ 
  int newstat=0;
  int status = DRMS_SUCCESS;
  DRMS_Record_t *outrec = NULL;
  DRMS_Segment_t *segout = NULL;
  DRMS_Array_t *outarr = NULL;
  DRMS_Type_t usetype = DRMS_TYPE_FLOAT;
  DRMS_RecLifetime_t lifetime;
  long long histrecnum=-1;
  int length[2];
  TIME trec, tnow, UNIX_epoch = -220924792.000; /* 1970.01.01_00:00:00_UTC */
  int ipsf=0;
  struct fresize_struct fress;
  static double defaultcoefs[] = {1.0, 0.459224, 0.132395, 0.019601, 0.000802, -4.31934E-05 };

  int errbufstat=setvbuf(stderr, NULL, _IONBF, BUFSIZ);
  int outbufstat=setvbuf(stdout, NULL, _IONBF, BUFSIZ);

/*
  double wt0, wt1, wt2, wt3, wt;
  double ut0, ut1, ut2, ut3, ut;
  double st0, st1, st2, st3, st;
  double ct0, ct1, ct2, ct3, ct;

  wt0=getwalltime();
  ct0=getcputime(&ut0, &st0);
*/

  char *outseries = (char *)cmdparams_save_str(&cmdparams, "out", &newstat);
  int verbflag = cmdparams_save_int(&cmdparams, "VERB", &newstat);
  int permflag = cmdparams_save_int(&cmdparams, "PERM", &newstat);
  if (permflag)
    lifetime = DRMS_PERMANENT;
  else
    lifetime = DRMS_TRANSIENT;

  char *histlinkname = (char *)cmdparams_save_str(&cmdparams, "histlink", &newstat);

  char *modefile = (char *)cmdparams_save_str(&cmdparams, "MODELIST", &newstat);
  int xpixels = cmdparams_save_int(&cmdparams, "XPIXELS", &newstat);
  int ypixels = cmdparams_save_int(&cmdparams, "YPIXELS", &newstat);
  double solarp = cmdparams_save_float(&cmdparams, "SOLARP", &newstat);
  double obsb0 = cmdparams_save_float(&cmdparams, "OBSB0", &newstat);
  double obsdist = cmdparams_save_float(&cmdparams, "OBSDIST", &newstat);
  float xoffset = cmdparams_save_float(&cmdparams, "XOFFSET", &newstat);
  float yoffset = cmdparams_save_float(&cmdparams, "YOFFSET", &newstat);
  int iinten = cmdparams_save_int(&cmdparams, "IINTEN", &newstat);
  float phasein = cmdparams_save_float(&cmdparams, "PHASE", &newstat);
  float freqin = cmdparams_save_float(&cmdparams, "FREQUENCY", &newstat);
  int ilimbdark = cmdparams_save_int(&cmdparams, "ILIMBDARK", &newstat);
  double coefs[6];
  int n_user_coefs = cmdparams_get_int(&cmdparams, "coefs_nvals", &status);
  if (ilimbdark)
  {
    if (n_user_coefs == 6)
    {
      double *cmdcoefs;
      int i;
      cmdparams_get_dblarr(&cmdparams, "coefs", &cmdcoefs, &status);
      for (i=0; i<6; i++)
        coefs[i] = cmdcoefs[i];
    }
    else
    {
      int i;
      for (i=0; i<6; i++)
        coefs[i] = defaultcoefs[i];
    }
  }
  int iheight = cmdparams_save_int(&cmdparams, "IHEIGHT", &newstat);

  int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
  int bxoff = cmdparams_save_int(&cmdparams, "BIN_XOFF", &newstat);
  int byoff = cmdparams_save_int(&cmdparams, "BIN_YOFF", &newstat);

  float sigma = cmdparams_save_float(&cmdparams, "GAUSS_SIG", &newstat);
  int hwidth = cmdparams_save_int(&cmdparams, "GAUSS_WID", &newstat);
  int nsub = cmdparams_save_int(&cmdparams, "GAUSS_NSUB", &newstat);
  int gxoff = cmdparams_save_int(&cmdparams, "GAUSS_XOFF", &newstat);
  int gyoff = cmdparams_save_int(&cmdparams, "GAUSS_YOFF", &newstat);
// for now use this patch so i can specify integers for GAUSS_SIG
  sigma/=sqrt(2);

/*
  int ibox = cmdparams_save_int(&cmdparams, "IBOX", &newstat);
  int igauss = cmdparams_save_int(&cmdparams, "IGAUSS", &newstat);
  float wpsfx = cmdparams_save_float(&cmdparams, "WPSFX", &newstat);
  float wpsfy = cmdparams_save_float(&cmdparams, "WPSFY", &newstat);
  int ilorentz = cmdparams_save_int(&cmdparams, "ILORENTZ", &newstat);
  float wlorentz = cmdparams_save_float(&cmdparams, "WLORENTZ", &newstat);
  int ilinint = cmdparams_save_int(&cmdparams, "ILININT", &newstat);
  int icubint = cmdparams_save_int(&cmdparams, "ICUBINT", &newstat);
  int ibin = cmdparams_save_int(&cmdparams, "IBIN", &newstat);
  int nbin = cmdparams_save_int(&cmdparams, "NBIN", &newstat);
  int iskip = cmdparams_save_int(&cmdparams, "ISKIP", &newstat);
  int noff = cmdparams_save_int(&cmdparams, "NOFF", &newstat);
*/
  double tstart = cmdparams_save_time(&cmdparams, "TSTART", &newstat);
  int dtoff = cmdparams_save_int(&cmdparams, "DTOFF", &newstat);
  int ndt = cmdparams_save_int(&cmdparams, "NDT", &newstat);

  double pangle = RADSINDEG*solarp;
  double b0 = RADSINDEG*obsb0;

//  if (ibox+igauss+ilorentz+ilinint+icubint > 0) ipsf=1;

  FILE *fptr = fopen(modefile,"r");
  if (fptr == NULL)
  {
    fprintf(stderr, "ERROR: can't open mode file %s\n",modefile);
    return 1;
  }

  if (newstat) 
  {
    fprintf(stderr, "ERROR: problem with input arguments, status = %d, diagnosis follows\n", newstat);
    cpsave_decode_error(newstat);
    return 1;
  }  
  else if (savestrlen != strlen(savestr)) 
  {
    fprintf(stderr, "ERROR: problem with savestr, savestrlen = %d, strlen(savestr) = %d\n", savestrlen, (int)strlen(savestr));
    return 1;
  }

  // set up ancillary dataseries for processing metadata
  DRMS_Record_t *tempoutrec = drms_create_record(drms_env, 
                                                 outseries,
                                                 DRMS_TRANSIENT, 
                                                 &status);

  if (status != DRMS_SUCCESS) 
  {
   fprintf(stderr,"ERROR: couldn't open a record in output dataseries %s, status = %d\n", outseries, status);
   return 1;
  }

  DRMS_Link_t *histlink = hcon_lookup_lower(&tempoutrec->links, histlinkname);
//  char *cvsinfo = strdup("$Header: /home/cvsuser/cvsroot/JSOC/proj/globalhs/apps/mkylms.c,v 1.5 2013/05/03 20:44:38 tplarson Exp $");
  if (histlink != NULL) 
  {
    histrecnum=set_history(histlink);
    if (histrecnum < 0)
    {
      drms_close_record(tempoutrec, DRMS_FREE_RECORD);
      return 1;
    }
  }
  else
  {
    fprintf(stderr,"WARNING: could not find history link in output dataseries\n");
  }


// these must be present in the output dataseries and variable, not links or constants   
  char *outchecklist[] = {"T_REC", "QUALITY", "T_OBS", "CRLT_OBS", "CRLN_OBS", 
                          //"SAT_ROT", "INST_ROT", "IM_SCALE",
                          "CDELT1", "CDELT2"};

  int itest;
  for (itest=0; itest < ARRLENGTH(outchecklist); itest++)
  {
    DRMS_Keyword_t *outkeytest = hcon_lookup_lower(&tempoutrec->keywords, outchecklist[itest]);
    if (!outkeytest || outkeytest->info->islink || outkeytest->info->recscope == 1)
    {
      fprintf(stderr, "ERROR: output keyword %s is either missing, constant, or a link\n", outchecklist[itest]);
      drms_close_record(tempoutrec, DRMS_FREE_RECORD);
      return 1;
    }
  }
  float cadence = drms_getkey_float(tempoutrec, "T_REC_step", &status);
  drms_close_record(tempoutrec, DRMS_FREE_RECORD);

  int i, hold;
  long nmodes=0;
  while ((hold = getc(fptr)) != EOF)
    if (hold == '\n')
      nmodes++;
  fclose(fptr);

  fptr = fopen(modefile,"r");
  int *lim, *mim;
  int lin, min;
  int mmax=0;
  int lmax=0;
  int *marr_lmax, *marr_count, **marr_index;
//  fscanf(fptr,"%d",&nmodes);
  lim=(int *)malloc(nmodes*sizeof(int));
  mim=(int *)malloc(nmodes*sizeof(int));
  for (i=0;i<nmodes;i++)
  {
    if (fscanf(fptr,"%d %d\n", &lin, &min) != 2)
    {
      fprintf(stderr,"ERROR: something went wrong on line i=%d\n", i);
      return 1;
    }
//    fscanf(fptr,"%d",&lin);
//    fscanf(fptr,"%d",&min);
    if (min > lin  || lin < 0)
    {
      fprintf(stderr,"ERROR: mode file has m > l or l < 0\n");
      return 1;
    }
/*
    if (min > mmax)
      mmax=min;
    if (lin > lmax)
      lmax=lin;
*/
    lim[i]=lin;
    mim[i]=min;
  }
  fclose(fptr);

  if (verbflag)
    printf("nmodes=%ld, dtoff=%d, ndt=%d\n", nmodes, dtoff, ndt); 

  int *lim2, *mim2;
  if (ndt > 0)
  {
    if (dtoff + ndt > nmodes)
      nmodes=nmodes-dtoff;
    else
      nmodes=ndt;

    lim2=(int *)malloc(nmodes*sizeof(int));
    mim2=(int *)malloc(nmodes*sizeof(int));
    for (i=0;i<nmodes;i++)
    {
      lim2[i]=lim[dtoff+i];
      mim2[i]=mim[dtoff+i];
    }
    free(lim);
    free(mim);
    lim=lim2;
    mim=mim2;
  }

  if (nmodes <= 0)
  {
    fprintf(stderr, "ERROR: no modes selected\n");
    return 1;
  }

  for (i=0;i<nmodes;i++)
  {
    if (mim[i] > mmax)
      mmax=mim[i];
    if (lim[i] > lmax)
      lmax=lim[i];
  }

  marr_lmax=(int *)(malloc((mmax+1)*sizeof(int)));
  marr_count=(int *)(malloc((mmax+1)*sizeof(int)));
  marr_index=(int **)(malloc((mmax+1)*sizeof(int *)));
  for (i=0;i<=mmax;i++)
  {
    marr_lmax[i]=-1;
    marr_count[i]=0;
    marr_index[i]=(int *)(malloc((lmax+1)*sizeof(int)));
  }
  for (i=0;i<nmodes;i++)
  {
    if (lim[i] > marr_lmax[mim[i]])
      marr_lmax[mim[i]]=lim[i];
    marr_index[mim[i]][marr_count[mim[i]]]=i;
    marr_count[mim[i]]++;
  }

//  double obsdist = 1.0;               //in au
  double distobs = obsdist*AU/RTRUE;  //in solar radii
  double dsunobs = obsdist*AU;        //in meters
  double x0=(double)xpixels/2-0.5 + xoffset;
  double y0=(double)ypixels/2-0.5 + yoffset;
  long npixels=xpixels*ypixels;
  double imscale=1.97784*1024/xpixels;
  double scale=imscale/(180*60*60/PI);
  length[0]=xpixels;
  length[1]=ypixels;
  float obsl0=211.67;
  int obscr=1784;
//  double rsunobs=atan(asin(RTRUE/AU/obsdist))*60*60/RADSINDEG;
//  double rsun=rsunobs/imscale;
  double rsunobs=asin(RTRUE/AU/obsdist)*60*60/RADSINDEG;
  double rsun=tan(asin(RTRUE/AU/obsdist))/scale;
  double cdelt=rsunobs/rsun;  // Use mean image scale across solar image
/*
; Sun is sitting at the center of the main coordinate system
; and has radius 1.
; Observer is at robs=(xobs,yobs,zobs) moving with a velocity
; vobs.
; Start by setting robs from distobs and b0
*/

  double robs_x = distobs * cos(b0);
  double robs_y = 0.0;
  double robs_z = distobs * sin(b0);

// Start with CCD coordinates
  double *x6 = (double *)(malloc(npixels*sizeof(double)));
  double *y6 = (double *)(malloc(npixels*sizeof(double)));
  int j;
  for (i=0;i<ypixels;i++)
    for (j=0;j<xpixels;j++)
      x6[i*ypixels+j]=(double)j;
  for (i=0;i<ypixels;i++)
    for (j=0;j<xpixels;j++)
      y6[i*ypixels+j]=(double)i;

  double *x4 = (double *)(malloc(npixels*sizeof(double)));
  double *y4 = (double *)(malloc(npixels*sizeof(double)));
  for (i=0;i<npixels;i++)
  {
    x4[i]=scale*(x6[i] - x0);
    y4[i]=scale*(y6[i] - y0);
  }


/*
; Rotate by the P-angle. New coordinate system has the y-axis pointing
; towards the solar north pole.
*/

  double *x2 = (double *)(malloc(npixels*sizeof(double)));
  double *y2 = (double *)(malloc(npixels*sizeof(double)));
  for (i=0;i<npixels;i++)
  {
    x2[i]= x4[i]*cos(pangle) + y4[i]*sin(pangle);
    y2[i]=-x4[i]*sin(pangle) + y4[i]*cos(pangle);
  }


/*
; Now transform to put the coordinates into the solar coordinate
; system.
; First find the directions (vecx and vecy) the x2 and y2 coordinate
; axis correspond to. vecsun points towards the Sun. Note that the
; (x2,y2,Sun) system is left handed. These vectors are unit vectors.
*/
  double vecx_x=0.0;
  double vecx_y=1.0;
  double vecx_z=0.0;
  double vecy_x=-sin(b0);
  double vecy_y=0.0;
  double vecy_z=cos(b0);
  double vecsun_x=-cos(b0);
  double vecsun_y=0.0;
  double vecsun_z=-sin(b0);
// Now the proper direction can be found. These are not unit vectors.
  double *x1 = (double *)(malloc(npixels*sizeof(double)));
  double *y1 = (double *)(malloc(npixels*sizeof(double)));
  double *z1 = (double *)(malloc(npixels*sizeof(double)));
  double *q1 = (double *)(malloc(npixels*sizeof(double)));
  for (i=0;i<npixels;i++)
  {
    x1[i]=vecx_x*x2[i] + vecy_x*y2[i] + vecsun_x;
    y1[i]=vecx_y*x2[i] + vecy_y*y2[i] + vecsun_y;
    z1[i]=vecx_z*x2[i] + vecy_z*y2[i] + vecsun_z;
    q1[i]=1/sqrt(x1[i]*x1[i] + y1[i]*y1[i] + z1[i]*z1[i]);
  }
// Make them unit vectors.
  for (i=0;i<npixels;i++)
  {
    x1[i]*=q1[i];
    y1[i]*=q1[i];
    z1[i]*=q1[i];
  }



// Now find intersection with the Sun.
// Solve quadratic equation |robs+q*[x1,y1,z1]|=1 for q
// a, b and c are terms in a*x^2+bx+c=0. a==1 since [x1,y1,z1] is unit vector.
//  double *b = (double *)(malloc(npixels*sizeof(double)));
//  double *d = (double *)(malloc(npixels*sizeof(double)));
  double b,d;
  double *q = (double *)(malloc(npixels*sizeof(double)));
  int *inflag = (int *)(malloc(npixels*sizeof(int)));
  double c = robs_x*robs_x + robs_y*robs_y + robs_z*robs_z -1;
  double *xsun = (double *)(malloc(npixels*sizeof(double)));
  double *ysun = (double *)(malloc(npixels*sizeof(double)));
  double *zsun = (double *)(malloc(npixels*sizeof(double)));
  double *cosrho = (double *)(malloc(npixels*sizeof(double)));
  for (i=0;i<npixels;i++)
  {
//    b[i]=2*(x1[i]*robs_x+y1[i]*robs_y+z1[i]*robs_z);
//    d[i]=b[i]*b[i] - 4*c;
    b=2*(x1[i]*robs_x+y1[i]*robs_y+z1[i]*robs_z);
    d=b*b - 4*c;
    if (d >= 0)
    {
      inflag[i]=1;
      q[i]=(-b - sqrt(d))/2;
      xsun[i]=robs_x + x1[i]*q[i];
      ysun[i]=robs_y + y1[i]*q[i];
      zsun[i]=robs_z + z1[i]*q[i];
      cosrho[i]=-(x1[i]*xsun[i]+y1[i]*ysun[i]+z1[i]*zsun[i]);
    }
    else
    {
      inflag[i]=0;
      xsun[i]=0;
      ysun[i]=0;
      zsun[i]=0;
      cosrho[i]=0;
    }
  }

  if (iheight)
  {
// Now calculate observing height in units of the solar radius
    double pc3[]={0.24047045,-0.76504325,0.86252178,-0.33859142};
    double height;
    double c0 = robs_x*robs_x + robs_y*robs_y + robs_z*robs_z;
    for (i=0;i<npixels;i++)
    {
      double x=1.0;
      height=pc3[0];
      for (j=1;j<4;j++)
      {
        x*=cosrho[i];
        height+=pc3[j]*x;
      }
      height*=1e6/RTRUE;
// Recalculate coordinates.
      b=2*(x1[i]*robs_x+y1[i]*robs_y+z1[i]*robs_z);
      c=c0-(1+height)*(1+height);
      d=b*b - 4*c;
      if (d >= 0)
      {
        inflag[i]=1;
        q[i]=(-b - sqrt(d))/2;
// Remember to divide by new target radius
        xsun[i]=(robs_x + x1[i]*q[i])/(1+height);
        ysun[i]=(robs_y + y1[i]*q[i])/(1+height);
        zsun[i]=(robs_z + z1[i]*q[i])/(1+height);
      }
      else
      {
        inflag[i]=0;
        xsun[i]=0;
        ysun[i]=0;
        zsun[i]=0;
      }
    }
  }

  double *phisun = (double *)(malloc(npixels*sizeof(double)));
  double *cphisun = (double *)(malloc(npixels*sizeof(double)));
  double *sphisun = (double *)(malloc(npixels*sizeof(double)));
  double *clatsun = (double *)(malloc(npixels*sizeof(double)));
  double *slatsun = (double *)(malloc(npixels*sizeof(double)));
  double *prad = (double *)(malloc(npixels*sizeof(double)));
  double *pphi = (double *)(malloc(npixels*sizeof(double)));
  double *plat = (double *)(malloc(npixels*sizeof(double)));

  for (i=0;i<npixels;i++)
  {
    phisun[i] = atan2(ysun[i],xsun[i]);
    cphisun[i] = cos(phisun[i]);
    sphisun[i] = sin(phisun[i]);
    slatsun[i] = zsun[i];
    clatsun[i] = sqrt(1 - zsun[i]*zsun[i]);
//  Set velocity projection factors
    prad[i]=xsun[i]*x1[i]+ysun[i]*y1[i]+zsun[i]*z1[i];
    pphi[i]=-sphisun[i]*x1[i]+cphisun[i]*y1[i];
    plat[i]=-cphisun[i]*slatsun[i]*x1[i]-sphisun[i]*slatsun[i]*y1[i]+clatsun[i]*z1[i];
  }

  double *ld;
  if (ilimbdark)
  {
    ld=(double *)malloc(npixels*sizeof(double));
    int iy,ix;
    for (iy=0; iy<ypixels; iy++)
      for (ix=0; ix<xpixels; ix++)
      {
        double costheta2;
        double xi, mu, z, ldi;
        double x, y, R2;
        int ord;
        int index=iy*xpixels + ix;

        if (!inflag[index])
          continue;

        /* get coordinates of point */
        x = ((double)ix - x0) / rsun; 
        y = ((double)iy - y0) / rsun;


        R2 = x*x + y*y;
        if (R2 <= 1.0)
          costheta2 = 1.0 - R2;
        else
          costheta2 = 0.0;

        mu = sqrt(costheta2);
	xi = log(mu);
	z = 1.0;
	ldi = 1.0;
	for (ord=1; ord<6; ord++)
        {
	  z *= xi;
	  ldi += coefs[ord] * z;
        }
        // not sure what to do with this...
        if (ldi <= 0.0 || isnan(ldi))
        {
          //*Op = missval;
          //ncropped++;
          ld[index]=0.0;
        }
	else
	  ld[index] = ldi;
	}
  }

  free(x6);free(y6);free(x4);free(y4);free(x2);free(y2);
  free(x1);free(y1);free(z1);free(q1);
//  free(b);free(d);
  free(cosrho);
  free(xsun);free(ysun);free(zsun);

  if (verbflag)
    printf("coordinate arrays set up\n");

  int lmax1, l, m, il;
  long nplm=10001;
  double *plm   = (double *)(malloc(nplm*(lmax+1)*sizeof(double)));
  double *dplm  = (double *)(malloc(nplm*(lmax+1)*sizeof(double)));
  double *plms  = (double *)(malloc(nplm*nmodes*sizeof(double)));
  double *dplms = (double *)(malloc(nplm*nmodes*sizeof(double)));
  double *dxplm = (double *)(malloc(nplm*sizeof(double)));
  for (i=0;i<nplm;i++)
    dxplm[i]= -1 + (double)i * 2 / (nplm - 1);
  int *indx = (int *)(malloc ((lmax+1)*sizeof(int)));
  for (l=0; l<=lmax; l++)
    indx[l]=l;

  for (m=0;m<=mmax;m++)
  {
    if (marr_lmax[m] == -1)
      continue;
    lmax1=marr_lmax[m];
    setplm2(m, lmax1, m, nplm, indx, dxplm, nplm, plm, dplm);
    for (i=0;i<marr_count[m];i++)
    {
      il=marr_index[m][i];
      l=lim[il];
      for (j=0;j<nplm;j++)
      {
        plms[il*nplm + j]=plm[l*nplm + j];
        dplms[il*nplm + j]=dplm[l*nplm + j];
      }
    }
  }

  free(plm);
  free(dplm);
  free(dxplm);

  if (verbflag)
    printf("plms set up\n");

//  float cadence=60.0;
//  double time0=612201600.0 + dtoff*cadence;
  double time0=tstart + dtoff*cadence;
  double *vradr=(double *)(malloc(npixels*sizeof(double)));
  double *vradi=(double *)(malloc(npixels*sizeof(double)));

//reuse arrays malloc'ed above to save memory
  double *vlatr, *vlati;
  vlatr=vradr;
  vlati=vradi;
//  double *vlatr=(double *)(malloc(npixels*sizeof(double)));
//  double *vlati=(double *)(malloc(npixels*sizeof(double)));
  double *vphir=(double *)(malloc(npixels*sizeof(double)));
  double *vphii=(double *)(malloc(npixels*sizeof(double)));

/*
  double *velr=(double *)(malloc(npixels*sizeof(double)));
  double *veli=(double *)(malloc(npixels*sizeof(double)));
  double *velsum=(double *)(malloc(npixels*sizeof(double)));
*/
  float *velr=(float *)(malloc(npixels*sizeof(float)));
  float *veli=(float *)(malloc(npixels*sizeof(float)));
  float *velsum=(float *)(malloc(npixels*sizeof(float)));
  float *velrsave=velr;
  float *velisave=veli;
  float *velssave=velsum;

  double *plminterp=(double *)(malloc(npixels*sizeof(double)));
  double ll, dinterp;

  float *binptr;
  int xpix1=xpixels;
  int ypix1=ypixels;
  if (nbin > 0)
  {
    x0=(x0 - bxoff + 0.5)/nbin - 0.5;
    y0=(y0 - byoff + 0.5)/nbin - 0.5;
//    imscale*=nbin;
    cdelt*=nbin;
    xpix1=xpixels/nbin;
    ypix1=ypixels/nbin;
    length[0]=xpix1;
    length[1]=ypix1;
    binptr=(float *)malloc(xpix1*ypix1*sizeof(float));
  }

  float *gaussptr;
  int xpix2, ypix2;
  if (hwidth > 0)
  {
    x0=(x0 - gxoff)/nsub;
    y0=(y0 - gyoff)/nsub;
//    imscale*=nsub;
    cdelt*=nsub;
    xpix2=xpix1/nsub;
    ypix2=ypix1/nsub;
    length[0]=xpix2;
    length[1]=ypix2;
    gaussptr=(float *)malloc(xpix2*ypix2*sizeof(float));  
    init_fresize_gaussian(&fress, sigma, hwidth, nsub);
  }

  double *phase=(double *)malloc(npixels*sizeof(double));
  double rsecs = RTRUE/C;
  for (i=0;i<npixels;i++)
  {
    if (inflag[i])
    {
     double deltat;
     deltat=(q[i]-distobs+1)*rsecs;
     phase[i]=phasein-2*PI*freqin*deltat;
    }
  }

  for (i=0;i<nmodes;i++)
  {
    trec=time0+i*cadence;
    if (verbflag > 1)
      printf("i=%d, trec=%f, l=%d, m=%d\n",i,trec,lim[i],mim[i]);

    outrec = drms_create_record(drms_env, outseries, lifetime, &status);
    if (status != DRMS_SUCCESS || outrec == NULL)
    {
      fprintf(stderr,"ERROR: unable to open record in output dataseries, i = %d, status = %d\n", i, status);
      return 0;
    }
    if (histlink)
      drms_setlink_static(outrec, histlinkname, histrecnum);

    for (j=0;j<npixels;j++)
      plminterp[j]=Ccintd(&plms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);

    for (j=0;j<npixels;j++)
    {
      if (inflag[j])
      {
        vradr[j]=1000 * (cos(mim[i]*phisun[j]+phase[j])) * plminterp[j];  //Ccintd(&plms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
        vradi[j]=1000 * (sin(mim[i]*phisun[j]+phase[j])) * plminterp[j];  //Ccintd(&plms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
// needed to make main leaks positive since DATASIGN=-1
        if (iinten)
        {
          if (ilimbdark)
          {
            velr[j]=-1*vradr[j]*ld[j];
            veli[j]=-1*vradi[j]*ld[j];
            velsum[j]=velr[j]+veli[j];
          }
          else
          {
            velr[j]=-1*vradr[j];
            veli[j]=-1*vradi[j];
            velsum[j]=velr[j]+veli[j];
          }
        }
        else
        {
          velr[j]=prad[j]*vradr[j];
          veli[j]=prad[j]*vradi[j];
          velsum[j]=velr[j]+veli[j];
        }
/*
if (j == 524800)
velsum[j]=1.0;
else
velsum[j]=0.0;
*/
      }
      else
      {
        velr[j]=0;
        veli[j]=0;
        velsum[j]=0;
      }
    }

    if (nbin > 0)
    {
      fbin(velr, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
      velr=binptr;
    }

    if (hwidth > 0)
    {
      fresize(&fress, velr, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
      velr=gaussptr;
    }

    segout = drms_segment_lookup(outrec, "vradr");
    outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velr, &status);
    if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
    {
      fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
      return 0;
    }
    drms_segment_write(segout, outarr, 0);
    free(outarr);

    if (nbin > 0)
    {
      fbin(veli, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
      veli=binptr;
    }

    if (hwidth > 0)
    {
      fresize(&fress, veli, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
      veli=gaussptr;
    }

    segout = drms_segment_lookup(outrec, "vradi");
    outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, veli, &status);
    if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
    {
      fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
      return 0;
    }
    drms_segment_write(segout, outarr, 0);
    free(outarr);

    if (nbin > 0)
    {
      fbin(velsum, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
      velsum=binptr;
    }

    if (hwidth > 0)
    {
      fresize(&fress, velsum, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
      velsum=gaussptr;
    }

    segout = drms_segment_lookup(outrec, "vradsum");
    outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velsum, &status);
    if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
    {
      fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
      return 0;
    }
    drms_segment_write(segout, outarr, 0);
    free(outarr);

    if (!iinten)
    {
      velr=velrsave;
      veli=velisave;
      velsum=velssave;
      for (j=0;j<npixels;j++)
      {
        if (inflag[j])
        {
          if (lim[i] != 0)
          {
            ll=sqrt(lim[i]*(lim[i]+1.0));
            dinterp=Ccintd(&dplms[i*nplm],nplm,(nplm-1)*(slatsun[j]+1)/2);
            vlatr[j]=1000 * (cos(mim[i]*phisun[j]+phase[j]))*dinterp/ll*clatsun[j];
            vlati[j]=1000 * (sin(mim[i]*phisun[j]+phase[j]))*dinterp/ll*clatsun[j];
            vphir[j]=1000 * mim[i]*(-sin(mim[i]*phisun[j]+phase[j]))*plminterp[j]/ll/clatsun[j];
            vphii[j]=1000 * mim[i]*( cos(mim[i]*phisun[j]+phase[j]))*plminterp[j]/ll/clatsun[j];
          }
          else
          {
            vlatr[j]=vlati[j]=0;
            vphir[j]=vphii[j]=0;
          }
          velr[j]=pphi[j]*vphir[j]+plat[j]*vlatr[j];
          veli[j]=pphi[j]*vphii[j]+plat[j]*vlati[j];
          velsum[j]=velr[j]+veli[j];
        }
        else
        {
          velr[j]=0;
          veli[j]=0;
          velsum[j]=0;
        }
      }

      if (nbin > 0)
      {
        fbin(velr, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
        velr=binptr;
      }

      if (hwidth > 0)
      {
        fresize(&fress, velr, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
        velr=gaussptr;
      }

      segout = drms_segment_lookup(outrec, "vhorr");
      outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velr, &status);
      if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
      {
        fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
        return 0;
      }
      drms_segment_write(segout, outarr, 0);
      free(outarr);

      if (nbin > 0)
      {
        fbin(veli, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
        veli=binptr;
      }

      if (hwidth > 0)
      {
        fresize(&fress, veli, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
        veli=gaussptr;
      }

      segout = drms_segment_lookup(outrec, "vhori");
      outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, veli, &status);
      if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
      {
        fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
        return 0;
      }
      drms_segment_write(segout, outarr, 0);
      free(outarr);

      if (nbin > 0)
      {
        fbin(velsum, binptr, xpixels, ypixels, xpixels, xpix1, ypix1, xpix1, nbin, bxoff, byoff, 0.0);
        velsum=binptr;
      }

      if (hwidth > 0)
      {
        fresize(&fress, velsum, gaussptr, xpix1, ypix1, xpix1, xpix2, ypix2, xpix2, gxoff, gyoff, 0.0);
        velsum=gaussptr;
      }

      segout = drms_segment_lookup(outrec, "vhorsum");
      outarr = drms_array_create(DRMS_TYPE_FLOAT, 2, length, velsum, &status);
      if (segout == NULL || outarr == NULL || status != DRMS_SUCCESS)
      {
        fprintf(stderr, "ERROR: problem with output segment or array, i = %d, status = %d\n", i, status);
        return 0;
      }
      drms_segment_write(segout, outarr, 0);
      free(outarr);

    }

    velr=velrsave;
    veli=velisave;
    velsum=velssave;

    drms_setkey_time(outrec, "T_REC", trec);
    drms_setkey_int(outrec, "L", lim[i]);
    drms_setkey_int(outrec, "M", mim[i]);
    drms_setkey_float(outrec, "XOFF", xoffset);
    drms_setkey_float(outrec, "YOFF", yoffset);
    drms_setkey_float(outrec, "PANGLE", solarp);
    drms_setkey_float(outrec, "BANGLE", obsb0);
//obviously these are redundant, but give series designer freedom to use either
    drms_setkey_float(outrec, "DISTOBS", distobs);
    drms_setkey_float(outrec, "OBSDIST", obsdist);

    drms_setkey_string(outrec, "FILE", modefile);

    drms_setkey_float(outrec, "CRPIX1", x0+1);
    drms_setkey_float(outrec, "CRVAL1", 0.0);
    drms_setkey_float(outrec, "CDELT1", cdelt);

    drms_setkey_float(outrec, "CRPIX2", y0+1);
    drms_setkey_float(outrec, "CRVAL2", 0.0);
    drms_setkey_float(outrec, "CROTA2", -solarp);
    drms_setkey_float(outrec, "CDELT2", cdelt);

    drms_setkey_time(outrec, "T_OBS", trec);
    drms_setkey_int(outrec, "QUALITY", 0);
    drms_setkey_float(outrec, "CADENCE", cadence);
    drms_setkey_float(outrec, "CRLT_OBS", obsb0);
    drms_setkey_float(outrec, "CRLN_OBS", obsl0);
    drms_setkey_int(outrec, "CAR_ROT", obscr);

    drms_setkey_double(outrec, "DSUN_OBS", dsunobs);
    drms_setkey_double(outrec, "RSUN_OBS", rsunobs);

    drms_setkey_double(outrec, "OBS_VR", 0.0);
    drms_setkey_double(outrec, "OBS_VW", 0.0);
    drms_setkey_double(outrec, "OBS_VN", 0.0);

    drms_setkey_int(outrec, "DATASIGN", -1);

    tnow = (double)time(NULL);
    tnow += UNIX_epoch;
    drms_setkey_time(outrec, "DATE", tnow);
    drms_close_record(outrec, DRMS_INSERT_RECORD);

  }

  if (hwidth > 0)
    free_fresize(&fress);

  if (verbflag)
    printf("module %s successful completion\n", cmdparams.argv[0]);

  return 0;
}



/* Calculate the interpolation kernel. */
void Ccker(double *u, double s)
{
   double s2, s3;
     
   s2= s * s;
   s3= s2 * s;
   u[0] = s2 - 0.5 * (s3 + s);
   u[1] = 1.5*s3 - 2.5*s2 + 1.0;
   u[2] = -1.5*s3 + 2.0*s2 + 0.5*s;
   u[3] = 0.5 * (s3 - s2);
}

/* Cubic convolution interpolation */
double Ccintd(double *f, int nx, double x)
{
   double  ux[4];
   /* Sum changed to double to speed up things */
   double sum;

   int ix, ix1, i;
     
   if (x < 1. || x >= (double)(nx-2))
     return 0.0;
     
   ix = (int)x;
   Ccker(ux,  x - (double)ix);
     
   ix1 = ix - 1;
   sum = 0.;
   for (i = 0; i < 4; i++)
      sum = sum + f[ix1 + i] * ux[i];
   return sum;
}