integr.c

Go to the documentation of this file.

 
#include "libtpcmodel.h"
/*****************************************************************************/
 
/*****************************************************************************/
int interpolate(
  double *x,
  double *y,
  int nr,
  double *newx,
  double *newy,
  double *newyi,
  double *newyii,
  int newnr
) {
  int i, j;
  double ty, tyi, tyii;
  double ox1, ox2, oy1, oy2, oi1, oi2, oii1, oii2, dt, ndt;
  int verbose=0;
 
  if(verbose>0) printf("in interpolate()\n");
  /* Check for data */
  if(nr<1 || newnr<1) return 1;
  /* Check that programmer understood that outp data must've been allocated */
  if(newy==NULL && newyi==NULL && newyii==NULL) return 2;
 
  /* Initiate first two input samples */
  ox1=ox2=x[0];
  oy1=oy2=y[0];
  /* Extrapolate the initial phase with triangle... */
  if(ox1>0.0) {
    oy1=0.0;
    /* unless:
       -first input sample y<=0
       -initial gap is longer than input TAC sampling frequency
    */
    if(y[0]>0.0 && (nr>1 && x[0]<=x[1]-x[0])) ox1=0.0;
  }
  /* Integrals too */
  oi1=oii1=0.0;
  oi2=oi1 + (ox2-ox1)*(oy1+oy2)/2.0;
  oii2=oii1 + (ox2-ox1)*(oi1+oi2)/2.0;
  if(verbose>2) {
    printf("ox1=%g oy1=%g oi1=%g oii1=%g\n", ox1, oy1, oi1, oii1);
    printf("ox2=%g oy2=%g oi2=%g oii2=%g\n", ox2, oy2, oi2, oii2);
  }
 
  /* Set interpolated data before input data (even imaginary) to zero */
  j=0; while(j<newnr && newx[j]<ox1) {
    ty=tyi=tyii=0.0; if(verbose>4) printf("  ndt=%g\n", ox1-newx[j]);
    if(verbose>4) printf("  j=%d newx=%g ty=%g tyi=%g tyii=%g\n", j, newx[j], ty, tyi, tyii);
    if(newy!=NULL) newy[j]=ty;
    if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii;
    j++;
  } 
 
  /* Set interpolated data between ox1 and ox2 */
  dt=ox2-ox1; if(dt>0.0) while(j<newnr && newx[j]<=ox2) {
    ndt=newx[j]-ox1; if(verbose>4) printf("  ndt=%g\n", ndt);
    ty=((oy2-oy1)/dt)*ndt + oy1; if(newy!=NULL) newy[j]=ty;
    tyi=oi1 + 0.5*(ty+oy1)*ndt; if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii= oii1 + 0.5*(tyi+oi1)*ndt;
    if(verbose>4) printf("  j=%d newx=%g ty=%g tyi=%g\n", j, newx[j], ty, tyi);
    j++;
  }
 
  /* Go through input data, sample-by-sample */
  for(i=0; i<nr && j<newnr; i++) {
 
    ox1=ox2; oy1=oy2; oi1=oi2; oii1=oii2;
    ox2=x[i]; oy2=y[i];
    oi2=oi1 + (ox2-ox1)*(oy1+oy2)/2.0;
    oii2=oii1 + (ox2-ox1)*(oi1+oi2)/2.0;
    if(verbose>>3) {
      printf("ox1=%g oy1=%g oi1=%g oii1=%g\n", ox1, oy1, oi1, oii1);
      printf("ox2=%g oy2=%g oi2=%g oii2=%g\n", ox2, oy2, oi2, oii2);
    }
    
    /* Calculate input sample distance */
    dt=ox2-ox1; if(dt<0.0) return 3; else if(dt==0.0) continue;
 
    /* Any need for interpolation between ox1 and ox2? */
    while(j<newnr && newx[j]<=ox2) {
      ndt=newx[j]-ox1;
      ty=((oy2-oy1)/dt)*ndt + oy1; if(newy!=NULL) newy[j]=ty;
      tyi=oi1 + 0.5*(ty+oy1)*ndt; if(newyi!=NULL) newyi[j]=tyi;
      if(newyii!=NULL) newyii[j]=tyii= oii1 + 0.5*(tyi+oi1)*ndt;
      if(verbose>5) printf("  j=%d newx=%g ty=%g tyi=%g\n", j, newx[j], ty, tyi);
      j++;
    }
 
  } // next input sample
 
  /* Set interpolated data after input data, assuming steady input */
  while(j<newnr) {
    ndt=newx[j]-ox2;
    ty=oy2;
    tyi=oi2 + oy2*ndt;
    tyii=oii2 + 0.5*(tyi+oi2)*ndt;
    if(newy!=NULL) newy[j]=ty;
    if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii;
    if(verbose>5) printf("  j=%d newx=%g ty=%g tyi=%g\n", j, newx[j], ty, tyi);
    j++;
  } 
 
  if(verbose>0) printf("out interpolate()\n");
  return 0;
}
 
int finterpolate(
  float *x,
  float *y,
  int nr,
  float *newx,
  float *newy,
  float *newyi,
  float *newyii,
  int newnr
) {
  int i, j;
  float ty, tyi, tyii;
  float ox1, ox2, oy1, oy2, oi1, oi2, oii1, oii2, dt, ndt;
  int verbose=0;
 
  if(verbose>0) printf("in finterpolate()\n");
  /* Check for data */
  if(nr<1 || newnr<1) return 1;
  /* Check that programmer understood that outp data must've been allocated */
  if(newy==NULL && newyi==NULL && newyii==NULL) return 2;
 
  /* Initiate first two input samples */
  ox1=ox2=x[0];
  oy1=oy2=y[0];
  /* Extrapolate the initial phase with triangle... */
  if(ox1>0.0) {
    oy1=0.0;
    /* unless:
       -first input sample y<=0
       -initial gap is longer than input TAC sampling frequency
    */
    if(y[0]>0.0 && (nr>1 && x[0]<=x[1]-x[0])) ox1=0.0;
  }
  /* Integrals too */
  oi1=oii1=0.0;
  oi2=oi1 + (ox2-ox1)*(oy1+oy2)/2.0;
  oii2=oii1 + (ox2-ox1)*(oi1+oi2)/2.0;
 
  /* Set interpolated data before input data (even imaginary) to zero */
  j=0; while(j<newnr && newx[j]<ox1) {
    ty=tyi=tyii=0.0;
    if(newy!=NULL) newy[j]=ty;
    if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii;
    j++;
  } 
 
  /* Set interpolated data between ox1 and ox2 */
  dt=ox2-ox1; if(dt>0.0) while(j<newnr && newx[j]<=ox2) {
    ndt=newx[j]-ox1;
    ty=((oy2-oy1)/dt)*ndt + oy1; if(newy!=NULL) newy[j]=ty;
    tyi=oi1 + 0.5*(ty+oy1)*ndt; if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii= oii1 + 0.5*(tyi+oi1)*ndt;
    j++;
  }
 
  /* Go through input data, sample-by-sample */
  for(i=0; i<nr && j<newnr; i++) {
 
    ox1=ox2; oy1=oy2; oi1=oi2; oii1=oii2;
    ox2=x[i]; oy2=y[i];
    oi2=oi1 + (ox2-ox1)*(oy1+oy2)/2.0;
    oii2=oii1 + (ox2-ox1)*(oi1+oi2)/2.0;
    
    /* Calculate input sample distance */
    dt=ox2-ox1; if(dt<0.0) return 3; else if(dt==0.0) continue;
 
    /* Any need for interpolation between ox1 and ox2? */
    while(j<newnr && newx[j]<=ox2) {
      ndt=newx[j]-ox1;
      ty=((oy2-oy1)/dt)*ndt + oy1; if(newy!=NULL) newy[j]=ty;
      tyi=oi1 + 0.5*(ty+oy1)*ndt; if(newyi!=NULL) newyi[j]=tyi;
      if(newyii!=NULL) newyii[j]=tyii= oii1 + 0.5*(tyi+oi1)*ndt;
      j++;
    }
 
  } // next input sample
 
  /* Set interpolated data after input data, assuming steady input */
  while(j<newnr) {
    ndt=newx[j]-ox2;
    ty=oy2;
    tyi=oi2 + oy2*ndt;
    tyii=oii2 + 0.5*(tyi+oi2)*ndt;
    if(newy!=NULL) newy[j]=ty;
    if(newyi!=NULL) newyi[j]=tyi;
    if(newyii!=NULL) newyii[j]=tyii;
    j++;
  } 
 
  if(verbose>0) printf("out finterpolate()\n");
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int integrate(
  double *x,
  double *y,
  int nr,
  double *yi
) {
  int j;
 
  if(nr==1 || x[0]<=(x[1]-x[0])) yi[0]=0.5*y[0]*x[0];
  else yi[0]=0; /*If the gap in the beginning is longer than time 
                 between first and second time point */
  for(j=1; j<nr; j++) yi[j]=yi[j-1]+0.5*(y[j]+y[j-1])*(x[j]-x[j-1]);
 
  return 0;
}
 
int fintegrate(
  float *x,
  float *y,
  int nr,
  float *yi
) {
  int j;
 
  if(nr==1 || x[0]<=(x[1]-x[0])) yi[0]=0.5*y[0]*x[0];
  else yi[0]=0; /*If the gap in the beginning is longer than time 
                 between first and second time point */
  for(j=1; j<nr; j++) yi[j]=yi[j-1]+0.5*(y[j]+y[j-1])*(x[j]-x[j-1]);
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int petintegrate(
  double *x1,
  double *x2,
  double *y,
  int nr,
  double *newyi,
  double *newyii
) {
  int i, allocated=0;
  double *ti, x, a;
 
 
  /* Check that data is increasing in time and is not (much) overlapping */
  if(nr<1 || x1[0]<0) return 1;
  for(i=0; i<nr; i++) if(x2[i]<x1[i]) return 2;
  for(i=1; i<nr; i++) if(x1[i]<=x1[i-1]) return 3;
 
  /* Allocate memory for temp data, if necessary */
  if(newyi==NULL) {
    allocated=1;
    ti=(double*)malloc(nr*sizeof(double)); if(ti==NULL) return 4;
  } else
    ti=newyi;
 
  /* Calculate the integral at the end of first frame */
  ti[0]=(x2[0]-x1[0])*y[0];
  /* If frame does not start at 0 time, add the area in the beginning */
  /* But only if the gap in the beginning is less than the lenght of the 
     first frame */
  if(x1[0]>0) {
    if(x1[0]<=x2[0]-x1[0]) {
       x=(x1[0]+x2[0])/2.0; 
       a=(x1[0]*(y[0]/x)*x1[0])/2.0; 
       ti[0]+=a;
    }
  }
 
  /* Calculate integrals at the ends of following frames */
  for(i=1; i<nr; i++) {
    /* Add the area of this frame to the previous integral */
    a=(x2[i]-x1[i])*y[i]; ti[i]=ti[i-1]+a;
    /* Check whether frames are contiguous */
    if(x1[i]==x2[i-1]) continue;
    /* When not, calculate the area of an imaginary frame */
    x=(x1[i]+x2[i-1])/2.0;
    a=(x1[i]-x2[i-1])*
      (y[i]-(y[i]-y[i-1])*(x2[i]+x1[i]-2.0*x)/(x2[i]+x1[i]-x2[i-1]-x1[i-1]));
    ti[i]+=a;
  }
 
  /* Copy integrals to output if required */
  if(allocated) for(i=0; i<nr; i++) newyi[i]=ti[i];
 
  /* Calculate 2nd integrals if required */
  if(newyii!=NULL) {
    newyii[0]=x2[0]*ti[0]/2.0;
    for(i=1; i<nr; i++)
      newyii[i]=newyii[i-1]+(x2[i]-x2[i-1])*(ti[i-1]+ti[i])/2.0;
  }
  
  /* Free memory */
  if(allocated) free((char*)ti);
 
  return 0;
}
 
int fpetintegrate(
  float *x1,
  float *x2,
  float *y,
  int nr,
  float *newyi,
  float *newyii
) {
  int i, allocated=0;
  float *ti, x, a;
 
 
  /* Check that data is increasing in time and is not (much) overlapping */
  if(nr<1 || x1[0]<0) return 1;
  for(i=0; i<nr; i++) if(x2[i]<x1[i]) return 2;
  for(i=1; i<nr; i++) if(x1[i]<=x1[i-1]) return 3;
 
  /* Allocate memory for temp data, if necessary */
  if(newyi==NULL) {
    allocated=1;
    ti=(float*)malloc(nr*sizeof(float)); if(ti==NULL) return 4;
  } else
    ti=newyi;
 
  /* Calculate the integral at the end of first frame */
  ti[0]=(x2[0]-x1[0])*y[0];
  /* If frame does not start at 0 time, add the area in the beginning */
  /* But only if the gap in the beginning is less than the lenght of the 
     first frame */
  if(x1[0]>0) {
    if(x1[0]<=x2[0]-x1[0]) {
       x=(x1[0]+x2[0])/2.0; 
       a=(x1[0]*(y[0]/x)*x1[0])/2.0; 
       ti[0]+=a;
    }
  }
 
  /* Calculate integrals at the ends of following frames */
  for(i=1; i<nr; i++) {
    /* Add the area of this frame to the previous integral */
    a=(x2[i]-x1[i])*y[i]; ti[i]=ti[i-1]+a;
    /* Check whether frames are contiguous */
    if(x1[i]==x2[i-1]) continue;
    /* When not, calculate the area of an imaginary frame */
    x=(x1[i]+x2[i-1])/2.0;
    a=(x1[i]-x2[i-1])*
      (y[i]-(y[i]-y[i-1])*(x2[i]+x1[i]-2.0*x)/(x2[i]+x1[i]-x2[i-1]-x1[i-1]));
    ti[i]+=a;
  }
 
  /* Copy integrals to output if required */
  if(allocated) for(i=0; i<nr; i++) newyi[i]=ti[i];
 
  /* Calculate 2nd integrals if required */
  if(newyii!=NULL) {
    newyii[0]=x2[0]*ti[0]/2.0;
    for(i=1; i<nr; i++)
      newyii[i]=newyii[i-1]+(x2[i]-x2[i-1])*(ti[i-1]+ti[i])/2.0;
  }
  
  /* Free memory */
  if(allocated) free((char*)ti);
 
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int interpolate4pet(
  double *x,
  double *y,
  int nr,
  double *newx1,
  double *newx2, 
  double *newy,
  double *newyi,
  double *newyii,
  int newnr
) {
  int ret, fi, overlap=0, zeroframe=0;
  double petx[3], pety[3], petyi[3], petyii[3], fdur;
  int verbose=0;
 
  if(verbose>0) printf("in interpolate4pet()\n");
 
  /* Check for data */
  if(nr<1 || newnr<1) return 1;
  /* Check that programmer understood that outp data must've been allocated */
  if(newy==NULL && newyi==NULL && newyii==NULL) return 2;
  /* Check that input data is not totally outside the input data */
  if(newx2[newnr-1]<=x[0] || newx1[0]>=x[nr-1]) return 3;
 
  /* Check frame lengths, also for overlap and zero length frames */
  for(fi=0; fi<newnr; fi++) {
    /* Calculate frame length */
    fdur=newx2[fi]-newx1[fi];
    /* If frame length is <0, that is an error */
    if(fdur<0.0) return 4;
    if(fdur==0.0) zeroframe++;
    /* Overlap? */
    if(fi>0 && newx2[fi-1]>newx1[fi]) overlap++;
  }
  if(verbose>1) {
    printf("overlap := %d\n", overlap);
    printf("zeroframe := %d\n", zeroframe);
  }
 
  /* Interpolate and integrate one frame at a time, if there is:
     -overlap in frames
     -frames of zero length
     -only few interpolated frames
     -if only integrals are needed
     -if neither of integrals is needed, then there is no temp memory to
      do otherwise
  */
  if(overlap>0 || zeroframe>0 || newnr<=3 || newy==NULL || (newyi==NULL && newyii==NULL) ) {
 
    if(verbose>1) printf("frame-by-frame interpolation/integration\n");
    for(fi=0; fi<newnr; fi++) {
      /* Set frame start, middle and end times */
      petx[0]=newx1[fi]; petx[2]=newx2[fi]; petx[1]=0.5*(petx[0]+petx[2]);
      /* Calculate frame length */
      fdur=petx[2]-petx[0];
      /* If frame length is <0, that is an error */
      if(fdur<0.0) return 4;
      /* If frame length is 0, then use direct interpolation */
      if(fdur==0.0) {
        ret=interpolate(x, y, nr, petx, pety, petyi, petyii, 1);
        if(ret) return 10+ret;
        if(newy!=NULL) newy[fi]=petyi[0];
        if(newyi!=NULL) newyi[fi]=petyi[0];
        if(newyii!=NULL) newyii[fi]=petyii[0];
        continue;
      }
      /* Calculate integrals at frame start, middle, and end */
      ret=interpolate(x, y, nr, petx, NULL, petyi, petyii, 3);
      if(ret) return 20+ret;
      /* Set output integrals, if required */
      if(newyi!=NULL) newyi[fi]=petyi[1];
      if(newyii!=NULL) newyii[fi]=petyii[1];
      /* Calculate frame mean, if required */
      if(newy!=NULL) newy[fi]=(petyi[2]-petyi[0])/fdur;
    } // next frame
 
  } else {
 
    if(verbose>1) printf("all-frames-at-once interpolation/integration\n");
    /* Set temp array */
    double *tp; if(newyii!=NULL) tp=newyii; else tp=newyi;
    /* Integrate at frame start times */
    ret=interpolate(x, y, nr, newx1, NULL, tp, NULL, newnr);
    if(ret) return 10+ret;
    /* Integrate at frame end times */
    ret=interpolate(x, y, nr, newx2, NULL, newy, NULL, newnr);
    if(ret) return 10+ret;
    /* Calculate average frame value */
    for(fi=0; fi<newnr; fi++)
      newy[fi]=(newy[fi]-tp[fi])/(newx2[fi]-newx1[fi]);
    /* Calculate integrals */
    if(newyi!=NULL || newyii!=NULL) {
      for(fi=0; fi<newnr; fi++) newx1[fi]+=0.5*(newx2[fi]-newx1[fi]);
      ret=interpolate(x, y, nr, newx1, NULL, newyi, newyii, newnr);
      if(ret) return 10+ret;
      for(fi=0; fi<newnr; fi++) newx1[fi]-=(newx2[fi]-newx1[fi]);
    }
  }
 
  if(verbose>0) printf("out interpolate4pet()\n");
  return 0;
}
 
int finterpolate4pet(
  float *x,
  float *y,
  int nr,
  float *newx1,
  float *newx2, 
  float *newy,
  float *newyi,
  float *newyii,
  int newnr
) {
  int ret, fi, overlap=0, zeroframe=0;
  float petx[3], pety[3], petyi[3], petyii[3], fdur;
  int verbose=0;
 
  if(verbose>0) printf("in finterpolate4pet()\n");
 
  /* Check for data */
  if(nr<1 || newnr<1) return 1;
  /* Check that programmer understood that outp data must've been allocated */
  if(newy==NULL && newyi==NULL && newyii==NULL) return 2;
  /* Check that input data is not totally outside the input data */
  if(newx2[newnr-1]<=x[0] || newx1[0]>=x[nr-1]) return 3;
 
  /* Check frame lengths, also for overlap and zero length frames */
  for(fi=0; fi<newnr; fi++) {
    /* Calculate frame length */
    fdur=newx2[fi]-newx1[fi];
    /* If frame length is <0, that is an error */
    if(fdur<0.0) return 4;
    if(fdur==0.0) zeroframe++;
    /* Overlap? */
    if(fi>0 && newx2[fi-1]>newx1[fi]) overlap++;
  }
  if(verbose>1) {
    printf("overlap := %d\n", overlap);
    printf("zeroframe := %d\n", zeroframe);
  }
 
  /* Interpolate and integrate one frame at a time, if there is:
     -overlap in frames
     -frames of zero length
     -only few interpolated frames
     -if only integrals are needed
     -if neither of integrals is needed, then there is no temp memory to
      do otherwise
  */
  if(overlap>0 || zeroframe>0 || newnr<=3 || newy==NULL || (newyi==NULL && newyii==NULL) ) {
 
    if(verbose>1) printf("frame-by-frame interpolation/integration\n");
    for(fi=0; fi<newnr; fi++) {
      /* Set frame start, middle and end times */
      petx[0]=newx1[fi]; petx[2]=newx2[fi]; petx[1]=0.5*(petx[0]+petx[2]);
      /* Calculate frame length */
      fdur=petx[2]-petx[0];
      /* If frame length is <0, that is an error */
      if(fdur<0.0) return 4;
      /* If frame length is 0, then use direct interpolation */
      if(fdur==0.0) {
        ret=finterpolate(x, y, nr, petx, pety, petyi, petyii, 1);
        if(ret) return 10+ret;
        if(newy!=NULL) newy[fi]=petyi[0];
        if(newyi!=NULL) newyi[fi]=petyi[0];
        if(newyii!=NULL) newyii[fi]=petyii[0];
        continue;
      }
      /* Calculate integrals at frame start, middle, and end */
      ret=finterpolate(x, y, nr, petx, NULL, petyi, petyii, 3);
      if(ret) return 20+ret;
      /* Set output integrals, if required */
      if(newyi!=NULL) newyi[fi]=petyi[1];
      if(newyii!=NULL) newyii[fi]=petyii[1];
      /* Calculate frame mean, if required */
      if(newy!=NULL) newy[fi]=(petyi[2]-petyi[0])/fdur;
    } // next frame
 
  } else {
 
    if(verbose>1) printf("all-frames-at-once interpolation/integration\n");
    /* Set temp array */
    float *tp; if(newyii!=NULL) tp=newyii; else tp=newyi;
    /* Integrate at frame start times */
    ret=finterpolate(x, y, nr, newx1, NULL, tp, NULL, newnr);
    if(ret) return 10+ret;
    /* Integrate at frame end times */
    ret=finterpolate(x, y, nr, newx2, NULL, newy, NULL, newnr);
    if(ret) return 10+ret;
    /* Calculate average frame value */
    for(fi=0; fi<newnr; fi++)
      newy[fi]=(newy[fi]-tp[fi])/(newx2[fi]-newx1[fi]);
    /* Calculate integrals */
    if(newyi!=NULL || newyii!=NULL) {
      for(fi=0; fi<newnr; fi++) newx1[fi]+=0.5*(newx2[fi]-newx1[fi]);
      ret=finterpolate(x, y, nr, newx1, NULL, newyi, newyii, newnr);
      if(ret) return 10+ret;
      for(fi=0; fi<newnr; fi++) newx1[fi]-=(newx2[fi]-newx1[fi]);
    }
 
  }
 
  if(verbose>0) printf("out finterpolate4pet()\n");
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int petintegral(
  double *x1,
  double *x2,
  double *y,
  int nr,
  double *ie,
  double *iie
) {
  int i;
  double x, last_x, last_x2, last_y, last_integral, box_integral, half_integral;
  double gap_integral, integral, integral2, frame_len, xdist, s;
 
  /* Check for data */
  if(nr<1 || nr<1) return(1);
  /* Check that programmer understood that output must've been allocated */
  if(ie==NULL && iie==NULL) return(2);
 
  /* Initiate values to zero */
  last_x=last_x2=last_y=last_integral=0.0;
  box_integral=integral=integral2=frame_len=0.0;
 
  for(i=0; i<nr; i++) {
    frame_len=x2[i]-x1[i]; if(frame_len<0.0) return(5);
    x=0.5*(x1[i]+x2[i]); xdist=x-last_x; 
    if(last_x>0.0 && xdist<=0.0) return(6);
    if(x<0) {
      if(ie!=NULL) ie[i]=integral;
      if(iie!=NULL) iie[i]=integral2;
      continue;
    }
    s=(y[i]-last_y)/xdist; /* slope */
    /*If there is a big gap in the beginning it is eliminated */
    if(i==0)if(x1[0]>x2[0]-x1[0]){last_x2=last_x=x1[0];}
    /*Integral of a possible gap between frames*/
    gap_integral=(x1[i]-last_x2)*(last_y+s*((last_x2+x1[i])/2.0-last_x)); 
    /*Integral from the beginning of the frame to the middle of the frame */
    half_integral=(x-x1[i])*(last_y+s*((x1[i]+x)/2.0-last_x));
    integral=box_integral+gap_integral+half_integral; 
    /* half_integral is not added to box because it is more accurate to 
       increment integrals of full frames */
    box_integral+=gap_integral+frame_len*y[i];
    integral2+=xdist*(integral+last_integral)*0.5;
    if(ie!=NULL) ie[i]=integral;
    if(iie!=NULL) iie[i]=integral2;
    last_x=x; last_x2=x2[i];
    last_y=y[i]; last_integral=integral;
  }
 
  return(0);
}
 
int fpetintegral(
  float *x1,
  float *x2,
  float *y,
  int nr,
  float *ie,
  float *iie
) {
  int i;
  float x, last_x, last_x2, last_y, last_integral, box_integral, half_integral;
  float gap_integral, integral, integral2, frame_len, xdist, s;
 
  /* Check for data */
  if(nr<1 || nr<1) return(1);
  /* Check that programmer understood that output must've been allocated */
  if(ie==NULL && iie==NULL) return(2);
 
  /* Initiate values to zero */
  last_x=last_x2=last_y=last_integral=0.0;
  box_integral=gap_integral=integral=integral2=frame_len=0.0;
 
  for(i=0; i<nr; i++) {
    frame_len=x2[i]-x1[i]; if(frame_len<0.0) return(5);
    x=0.5*(x1[i]+x2[i]); xdist=x-last_x; if(last_x>0.0 && xdist<=0.0) return(6);
    if(x<0) {
      if(ie!=NULL) ie[i]=integral;
      if(iie!=NULL) iie[i]=integral2;
      continue;
    }
    s=(y[i]-last_y)/xdist; /* slope */
    /*If there is a big gap in the beginning it is eliminated */
    if(i==0)if(x1[0]>x2[0]-x1[0]){last_x2=last_x=x1[0];}
    /*Integral of a possible gap between frames*/
    gap_integral=(x1[i]-last_x2)*(last_y+s*((last_x2+x1[i])/2.0-last_x));
    /*Integral from the beginning of the frame i to the middle of the frame */
    half_integral=(x-x1[i])*(last_y+s*((x1[i]+x)/2.0-last_x));
    integral=box_integral+gap_integral+half_integral; 
    /* half_integral is not added to box because it is more accurate to 
       sum integrals of full frames */
    box_integral+=gap_integral+frame_len*y[i];
    integral2+=xdist*(integral+last_integral)*0.5;
    if(ie!=NULL) ie[i]=integral;
    if(iie!=NULL) iie[i]=integral2;
    last_x=x; last_x2=x2[i];
    last_y=y[i]; last_integral=integral;
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int petintegrate2fe(
  double *x1,
  double *x2,
  double *y,
  int nr,
  double *e,
  double *ie,
  double *iie
) {
  int i;
  double x, last_x, last_x2, last_y, last_integral;
  double value, integral, integral2, frame_len, xdist, s;
 
  /* Check for data */
  if(nr<1 || nr<1) return(1);
  /* Check that programmer understood that output must've been allocated */
  if(e==NULL && ie==NULL && iie==NULL) return(2);
 
  /* Initiate values to zero */
  last_x=last_x2=last_y=last_integral=value=integral=integral2=frame_len=s=0.0;
 
  for(i=0; i<nr; i++) {
    frame_len=x2[i]-x1[i]; if(frame_len<0.0) return(5);
    x=0.5*(x1[i]+x2[i]); xdist=x-last_x; if(last_x>0.0 && xdist<=0.0) return(6);
    if(x<0) {
      if(e!=NULL) e[i]=value;
      if(ie!=NULL) ie[i]=integral;
      if(iie!=NULL) iie[i]=integral2;
      continue;
    }
    s=(y[i]-last_y)/xdist; /* slope between x[i-1] and x[i] */
    /* If there is a big gap in the beginning, it is eliminated */
    if(i==0 && x1[0]>x2[0]-x1[0]) {last_x2=last_x=x1[0];}
    integral+=(x1[i]-last_x2)*(last_y+s*((last_x2+x1[i])/2.0-last_x)); /*gap*/
    integral+=frame_len*y[i];
    integral2+=(x2[i]-last_x2)*(integral+last_integral)*0.5;
    if(e!=NULL && i>0) {
      value=last_y+s*(last_x2-last_x); /* value at previous frame end */
      e[i-1]=value;
    }
    if(ie!=NULL) ie[i]=integral;
    if(iie!=NULL) iie[i]=integral2;
    last_x=x; last_x2=x2[i]; last_y=y[i]; last_integral=integral;
  }
  if(e!=NULL) {
    value=last_y+s*(last_x2-last_x); /* Value for the last frame */
    e[i-1]=value;
  }
 
  return(0);
}
 
int fpetintegrate2fe(
  float *x1,
  float *x2,
  float *y,
  int nr,
  float *e,
  float *ie,
  float *iie
) {
  int i;
  float x, last_x, last_x2, last_y, last_integral;
  float value, integral, integral2, frame_len, xdist, s;
 
  /* Check for data */
  if(nr<1 || nr<1) return(1);
  /* Check that programmer understood that output must've been allocated */
  if(e==NULL && ie==NULL && iie==NULL) return(2);
 
  /* Initiate values to zero */
  last_x=last_x2=last_y=last_integral=value=integral=integral2=frame_len=s=0.0;
 
  for(i=0; i<nr; i++) {
    frame_len=x2[i]-x1[i]; if(frame_len<0.0) return(5);
    x=0.5*(x1[i]+x2[i]); xdist=x-last_x; if(last_x>0.0 && xdist<=0.0) return(6);
    if(x<0) {
      if(e!=NULL) e[i]=value;
      if(ie!=NULL) ie[i]=integral;
      if(iie!=NULL) iie[i]=integral2;
      continue;
    }
    s=(y[i]-last_y)/xdist; /* slope between x[i-1] and x[i] */
    /* If there is a big gap in the beginning, it is eliminated */
    if(i==0 && x1[0]>x2[0]-x1[0]) {last_x2=last_x=x1[0];}
    integral+=(x1[i]-last_x2)*(last_y+s*((last_x2+x1[i])/2.0-last_x)); /*gap*/
    integral+=frame_len*y[i];
    integral2+=(x2[i]-last_x2)*(integral+last_integral)*0.5;
    if(e!=NULL && i>0) {
      value=last_y+s*(last_x2-last_x); /* value at previous frame end */
      e[i-1]=value;
    }
    if(ie!=NULL) ie[i]=integral;
    if(iie!=NULL) iie[i]=integral2;
    last_x=x; last_x2=x2[i]; last_y=y[i]; last_integral=integral;
  }
  if(e!=NULL) {
    value=last_y+s*(last_x2-last_x); /* Value for the last frame */
    e[i-1]=value;
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/