#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include "myc_fortran.h"

void F_FUN(dat_axis_range)(float *xin, int *nin, float *xstart, float *xend, float *pMargin, int *i1, int *i2) {
  int i;
  float x1, x2, margin;

  margin=*pMargin;
  if (margin<0) margin=-margin;
  if (*xstart>*xend) {
    x1=*xend-margin;
    x2=*xstart+margin;
  } else {
    x1=*xstart-margin;
    x2=*xend+margin;
  }
  *i1=*nin;
  *i2=0;
  for (i=0; i<*nin; i++) {
    if (x1 < xin[i] && xin[i] < x2) {
      if (i>*i2) *i2=i;
      if (i<*i1) *i1=i;
    }
  }
}

void dat_fillin(float *src, float *dst, int size, float *wsum, int j, float weight) {
  int k;

  dst += j*size;
  if (wsum!=NULL) wsum[j] += weight;
  for (k=0; k<size; k++) {
    *dst += weight * *src;
    dst++; src++;
  }
}

void dat_compress(float *xin, float *src, int nin, int grpsize, int nblocks
  , float xstart, float xs, int nout, float *dst, float *wout) {
/*
** compress (sum up) nin input channels to nout output channels
** source and destination may be multidimensional arrays, but the compression
** is performed only in one dimension (d). The shape of the array is summarized
** by the following values:
**   grpsize  product of the size of all dimensions lower (d)
**   nblocks  product of the size of all dimensions higher (d)
**   nin      length of compressed dimension
**    Examples:
**       number of dimensions        1       2            3
**       x-range                  0..9    0..9         0..4
**       y-range                          0..2         0..6
**       z-range                                       0..8
**       index formula               x  x+y*10  x+y*(5+z*7)
**       dimension to compress       x       x            y
**       grpsize                     1       1            5
**       nblocks                     1       3            9
**       nin                        10      10            7
**
** geometric principe of compression:
**       x[0]        x[1]        x[2]    x[3]  x[4]  x[5]
**        |           |           |       |     |     |
**  |-----0-----|-----1-----|-----2---|---3--|--4--|--5--|
**            |----0----|----1----|----2----|----3----|
**                 |         |         |         |     
**              x-start      |<--xs--->|
**
** The values from the source are distributed to the destination according to the overlap
** of the intervals. Output channel 0 will get a small part from input channel 0 and
** a large part from channel 1, output channel 5 will get all of channel 4 plus
** some part of channel 3 and channel 5.
**
** xin      x-values (size nin) (may be NULL, if x-values are [0.0,1.0,2.0,3.0,4.0,....])
** src      source (size grpsize*nin*nblocks)
** xstart   start value
** xs       x-step for output
** nout     number of output channels
** dst      destination (size grpsize*nout*nblocks)
** wout     output weights (size nout, may be NULL if not needed)
*/
  int i, j, m, ja, jb, bsize;
  float w, x, xl, xr, xa, xb, dx0, dx1, dx;

  assert(nin>=1 && src !=NULL && grpsize>0 && nblocks>0 && xs!=0.0 && nout>0 && dst!=NULL && wout!=NULL);
  for (j=0; j<nout; j++) {
    wout[j]=0;
  }
  bsize=grpsize*nout;
  for (j=0; j<bsize*nblocks; j++) {
    dst[j]=0;
  }
  for (m=0; m<nblocks; m++) {
    if (xin==NULL) {
      dx=1.0/xs;
      dx0=dx/2;
    } else {
      dx0=(xin[1]-xin[0])/xs*0.5;
    }
    if (nout==1) {
      for (i=0; i<nin; i++) {
        if (xin==NULL) {
          xa=(i-xstart)/xs-dx0;
        } else {
          if (i<nin-1) dx1=(xin[i+1]-xin[i])/xs*0.5;
          dx=dx0+dx1;
          xa=(xin[i]-xstart)/xs-dx0;
        }
        xb=xa+dx;
        if (xa>-0.5) {
          if (xb<0.5) {
            w=1.0;
          } else {
            w=(0.5-xa)/dx;
          }
        } else {
          if (xb<0.5) {
            w=(xb+0.5)/dx;
          } else {
            w=1.0;
          }
        }
        if (w>0.0) {
          dat_fillin(src, dst, grpsize, wout, 0, w);
        }
        src+=grpsize;
      }
    } else {
      for (i=0; i<nin; i++) {
        if (xin==NULL) {
          xa=(i-xstart)/xs-dx0;
        } else {
          if (i<nin-1) dx1=(xin[i+1]-xin[i])/xs*0.5;
          dx=dx0+dx1;
          xa=(xin[i]-xstart)/xs-dx0;
        }
        ja=xa; if (ja>xa) ja--;
        xb=xa+dx;
        jb=xb; if (jb<xb) jb++;
        if (ja<0) {
          ja=-1;
          if (xa<-1.0) xa=-1.0;
        }
        if (jb>=nout) {
          jb=nout;
          if (xb>nout) xb=nout;
        }   
        if (ja<jb) {
          if (ja+1==jb) {
            x=(xa+xb)*0.5-ja;
            w=(xb-xa)/dx;
            if (ja>=0) dat_fillin(src, dst, grpsize, wout, ja, (1-x)*w);
            if (jb<nout) dat_fillin(src, dst, grpsize, wout, jb, x*w);
          } else {
            xl=ja+1-xa;
            xr=xb+1-jb;
            w=1.0/dx;
            if (ja>=0) dat_fillin(src, dst, grpsize, wout, ja, xl*xl*w/2);
            if (ja+2==jb) {
              dat_fillin(src, dst, grpsize, wout, ja+1, ((1-xl/2)*xl+(1-xr/2)*xr)*w);
            } else {
              dat_fillin(src, dst, grpsize, wout, ja+1, ((1-xl/2)*xl+0.5)*w);
              for (j=ja+1; j<jb-1; j++) {
                dat_fillin(src, dst, grpsize, wout, j, w);
              }
              dat_fillin(src, dst, grpsize, wout, jb-1, ((1-xr/2)*xr+0.5)*w);
            }
            if (jb<nout) dat_fillin(src, dst, grpsize, wout, jb, xr*xr*w/2);
          }
        }
        src+=grpsize;
      }
    }
    wout=NULL; /* wout is stored only in the first block */
    dst+=bsize;
  }
}

int F_FUN(dat_comp)(int *flgs, float *xin, float **fin, int *nin, int *grpsize, int *nblocks
  , float *xstart, float *xs, int *nout, float **weights) {
/* flgs 0,2: xin unused      1,3: xin used
** flgs 0,1: weights unused  2,3: weights used
*/
  float *fout, *wout;
  int l;

  wout=calloc(*nout, sizeof(float));
  if (wout==NULL) return 0;
  l = *nout * *grpsize * *nblocks;
  fout=calloc(l, sizeof(float));
  if (fout==NULL) { free(wout); return 0; }
  if (0==(*flgs & 1)) xin=NULL;
  dat_compress(xin, *fin, *nin, *grpsize, *nblocks, *xstart, *xs, *nout, fout, wout);
  free(*fin);
  *fin=fout;
  if (0==(*flgs &2) || weights==NULL) {
    free(wout);
  } else {
    *weights=wout;
  }
  return l;
}

void F_FUN(dat_copy2f)(float **fin, float *fout, int *size, int *freeIt) {
  float *p;
  int i, s;

  p=*fin;
  s=*size;
  for (i=0; i<s; i++) {
    *fout++ = *p++;
  }
  if (*freeIt) {
    free(*fin);
    *fin=NULL;
  }
}

void F_FUN(dat_copy2p)(float **fout, float *fin) {
  float *p;
  int i, s;

  p=calloc(1,sizeof(float));
  if (p==NULL) return;
  *p=*fin;
  *fout = p;
  return;
}

void F_FUN(dat_extract)(float **fin, int *offset, int *step, float *fout, int *size, int *freeIt) {
  float *p;
  int i, s;

  p=*fin + *offset;
  s=*size;
  for (i=0; i<s; i++) {
    *fout++ = *p;
    p += *step;
  }
  if (*freeIt) {
    free(*fin);
    *fin=NULL;
  }
}

#include "napi.h"

int F_FUN(dat_nexus_getslab)(NXhandle handle, int *pType, int *rank,
                             int *start, int *size, float **fout) {

  int status, i, type;
  int cellsize, ncells;
  void *data;
  float *flts;
  double *dptr;
  int *ip;
  int cstart[32], csize[32];

  type=*pType;
  switch (type) {
  case NX_FLOAT64: cellsize=8; break;
  case NX_UINT32: type=NX_INT32; /* fall through */
  case NX_INT32: /* fall through */
  case NX_FLOAT32: cellsize=4; break;
  case NX_UINT16: type=NX_INT16; /* fall through */
  case NX_INT16: cellsize=2; break;
  case NX_UINT8: type=NX_INT8; /* fall through */
  case NX_INT8: /* fall through */
  case NX_CHAR: cellsize=1; break;
  }
  if (cellsize<4 || *rank<1 || *rank>32) return NX_ERROR; /* for now, support only f32, i32 and f64 */
  ncells=1;
  for (i=0; i<*rank; i++) {
    cstart[i]=start[*rank-i-1];
    csize[i]=size[*rank-i-1];
    ncells=ncells*size[i];
  }
  for (;i<32;i++) {
    cstart[i]=0;
    csize[i]=1;
  }
  data=calloc(ncells, cellsize);
  if (data==NULL) return NX_ERROR;
  flts=(float *)data;

  status=NXgetslab(handle, data, cstart, csize);
  if (status!=NX_OK) { free(data); return NX_ERROR; }

  if (type==NX_INT32) { /* convert int32 to float 32 */
    ip=(int *)data;
    for (i=0; i<ncells; i++) {
      flts[i]=*ip; ip++;
    }
  } else if (type==NX_FLOAT64) { /* convert float 64 to float 32 */
    dptr=(double *)data;
    for (i=0; i<ncells; i++) {
      flts[i]=*dptr; dptr++;
    }
  }
  *fout=flts;
  return NX_OK;
}