/* Program written by Mark Keir to emulate satplt */
/* Language used: C */
/* Date of last revision: 20-02-97 */

/* Parameter definitions */

#define MAXN 92400   /* maximum number of points in data array */
#define MAXL 31      /* maximum length for cumulative sum */
#define MINW 20      /* minimum mark length */
#define MAXW 3000    /* maximum mark length */
#define AVLEN 21     /* average mark length */
#define NSUSP 2      /* number of suspect points */
#define LMAX 21      /* length for mean and std */
#define T 3.0        /* t ratio */
#define TOUT 20.0    /* t ratio for outliers */
#define DOUT 5.0     /* minimum threshold for outliers */ 
#define SDMAX 8.0    /* maximum allowable std in mark */
#define SDMIN 12.0   /* minimum possible difference for mark */
#define MAXREJ 4     /* maximum number of over-undershoots */
#define MAXOUT 3     /* maximum number of consecutive outliers in mark */
#define DMAX 5.0     /* threshold for over-undershoots */
#define MAXM 350     /* moving average filter length */

/* Libraries includes for functions */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Global variable declarations */

double y[MAXN];
double w[MAXN];
int im1,im2;
int nx;
double xs;
double xss;
int q;
int isv; /* 0:normal 1:possible mark 2:mark 3:possible normal -1:unknown*/
int nrej; /* number of rejected points in transition */
int nout; /* number of outliers in mark */
double xp;
double xb; /* previous point, point before possible mark */

/* User defined functions */

FILE *file_open(char name[],char access_mode[]);
void rmmarks(int n);
void adjmrk(int n);
int nxmark(double x);
void cumul(int lmax,double x);
double dmean(void);
double stdev(void);
double amean(int ns,int nf);
double yinter(double x1,double x2,double y1,double y2,double xi);
void plotcntrl(char head[]);
void filtdath(int k);
void filtdatl(int k);
/* line 51 */
/* Main program body */

main(int argc,char *argv[])   
{
/* File pointers */

 FILE *fin;
 FILE *fout;
 
/* Variable definitions */ 

 char filename[13];
 char header[87];
 char opt[8],mangle,seg;
 int i,j,k;
 int nrs,nrf;
 int ns[2];
 float tt;

 /* Actual file processing */
 
 strcpy(filename,argv[1]);
 fin=file_open(filename,"rb");
 
 /* Loop to read header */
 
 for (j=1;j<87;j++)
 {header[j]=getc(fin);
  if(header[j]<33)
  header[j]=' ';}
 
 /* Loop to read in data */
 
 i=1;
 while(((y[i]=getc(fin))!=EOF)&&(i<=MAXN))
 {i++;}
 fclose(fin);
 
 /* Read input variables */
 /*printf("Enter start time and finish time (in seconds):\n");*/
 /*scanf("%d %d",&ns[1],&ns[2]);*/
 /*printf("Time mark removal option (R,F,N):\n");*/
 /*scanf("%c",&mangle);*/
 sscanf(argv[2],"%c",&seg);
 sscanf(argv[3],"%c",&mangle);
 
 switch(seg){
 case 'S':
   /*start of record*/
 ns[1]=0;ns[2]=400;
 break;
 case 'M':
   /*middle of record*/
 ns[1]=400;ns[2]=800;
 break;
 case 'E':
   /*end of record*/
 ns[1]=800;ns[2]=1320;
 break;
 case 'A':
   /*all of the record*/
 ns[1]=0;ns[2]=1320;
 break;}
 
 nrs=ns[1]*70;
 nrf=ns[2]*70;
 /* printf("nrf=%d\n",nrf);*/
 if(nrs<1)
 {nrs=1;}
 if(nrf>MAXN)
 {nrf=MAXN;}

 switch(mangle){
  case 'R':
  /* marks removed only */
  rmmarks(i);
  break;
  case 'F':
  /* marks removed and filtered */
  rmmarks(i);
  filtdath(i);
  break;
  case 'N':
  /* raw data */
  break;}
 
 /* Output data to file for examination or plotting */

 plotcntrl(header);
 /* printf("nrf=%d\n",nrf);*/
 fout=fopen("sat_plot.dat","w");
 for(k=nrs;k<nrf;k++) 
 {tt=k/70.0;
  if((k%5)==1)
  {fprintf(fout," %g %g\n",tt,y[k]);}}
 fclose(fout);
 /* End program */
 
 return 0;
}

/* Function to verify data file existence */

FILE *file_open(char name[],char access_mode[])
{
 FILE *f;
 f=fopen(name,access_mode);
 if (f==NULL)
 {perror("Cannot open file");
  exit(EXIT_FAILURE);}
 return f;
}

/* Function to remove time marks from data stream */

void rmmarks(int n)  
{
 int c;
 double xv;
/* Initialise variables */
 for(c=1;c<n;c++)
 {
  if(c<=n)
  {xv=y[c];}
  else
  {xv=-1;}
  if(nxmark(xv)==1)
  {adjmrk(n);}
 }
 return;
}


int nxmark(double x)
{
 double xm,xsd;
 int nxtmrk;
 double maxv,maxv2;
 
 nxtmrk=0;
 if(x==-1.0)
 {
  if(isv==2)
  {im2=q;}
  if((isv==2)||(isv==3))
  {nxtmrk=1;}
  else
  {nxtmrk=0;}
  return nxtmrk;
 }
 if(im1==0)
 {
  isv=-1;
  im1=1;
  q=0;
  nx=0;
 }
 q++;
 xm=dmean();
 xsd=stdev();
 if((isv==1)&&(nx==1)) /* possible beginning */
 {
  if(fabs(xm-x)<=SDMAX) /* assumed mark */
  {isv=2;}
  else    
  {
   if((nrej>MAXREJ)||((xb-x)<=SDMIN))     /* not mark */
   {isv=0;}
   else  /* possible transition */
   {nx=0;
    nrej++;
   }
  }
 }
 if((isv==3)&&(nx==1)) /* possible end */
 {
  if((nrej<=MAXREJ)&&(fabs(x-xm)>DMAX))  /* reject point */
  {nx=0;
   im2=q-1;
   nrej++;
  }
  else  /* accept as mark */
  {isv=0;
   nxtmrk=1;
  }
 }
 if((T*xsd)>SDMIN)
 {maxv=T*xsd;}
 else
 {maxv=SDMIN;}
 if(nx>0)
 {
  if(((isv==0)||(isv==-1))&&((xp-x)>SDMIN)&&((xm-x)>maxv)) /* possible mark beginning */
  {isv=1;
    nx=0;
    xb=xp;  /* save point before mark */
    nrej=0;
    im1=q;
    im2=0;}
  else if(((isv==2)||(isv==-1))&&(xsd>SDMAX)) /* not in mark- too much variation */
  {isv=0;}
  else if(((isv==2)||(isv==-1))&&((q-im1)>MAXW)) /* not in mark- too long */
  {isv=0;}
  else if(((isv==2)||(isv==-1))&&((x-xm)>maxv)) 
  {
   if(fabs(q-im1)<MINW) /* not in mark- too short */
     {isv=0;}
   else   /* possible point after end mark */
   {isv=3;
    nx=0;
    nrej=0;
    im2=q-1;
   }
  }
 }                          /* line 252 */
 if((TOUT*xsd)>DOUT)
 {maxv2=TOUT*xsd;}
 else
 {maxv2=DOUT;}
 if((isv==2)&&(fabs(x-xm)>=maxv2))   /* outlier */
  {nout++;
  if(nout>MAXOUT) /* shift in mark level */
  {isv=0;         /* rejects mark */
    nx=0;         /* starts a new series for mean and std */
    cumul(LMAX,x);
  }
 }
 else             /* normal point */
 {nout=0;
  xp=x;
  cumul(LMAX,x);
 } 
 /*printf("%g\n",fabs(x-xm));*/
return nxtmrk;
}

/* Function to change datastream values within time marks */

void adjmrk(int n)  
{
 double shift;
 int i,j,na,ib,ie;
 double v[3];

 
 if(im1-NSUSP-AVLEN>1)
 {ib=im1-NSUSP-AVLEN;}
 else
 {ib=1;}
 na=im1-NSUSP-ib;
 v[1]=amean(ib,(ib+na));
 if(((im1+im2)/2-(AVLEN-1)/2)>(im1+NSUSP+1))
 {ib=(im1+im2)/2-(AVLEN-1)/2;}
 else
 {ib=im1+NSUSP+1;}
 if((im2-ib-NSUSP)<AVLEN)
 {na=im2-ib-NSUSP;}
 else
 {na=AVLEN;}
 v[2]=amean(ib,(ib+na));
 ib=im2+NSUSP+1;
 if((n-ib+1)<AVLEN)
 {na=n-ib-1;}
 else
 {na=AVLEN;}        /* line 302 */
 v[3]=amean(ib,(ib+na));
 if((im1>1)&&(im2>=n))
 {shift=v[1]-v[2];}
 else if((im1<=1)&&(im2<n))
 {shift=v[3]-v[2];}
 else if((im1>1)&&(im2<n))
 {shift=(v[1]+v[3])/2.0-v[2];}
 else
 {im1=1;im2=n;}
 for(i=im1;i<=im2;i++)
 {y[i]=y[i]+shift;}
 if((im1-NSUSP)>2)
  {ib=im1-NSUSP;}
 else
  {ib=2;}
 if((im1+NSUSP)<(n-1))
  {ie=im1+NSUSP;}
 else
  {ie=n-1;}
 for(i=ib;i<=ie;i++)
  {y[i]=yinter((ib-1),(ie+1),y[(ib-1)],y[(ie+1)],i);}
 if((im2-NSUSP)>2)
  {ib=im2-NSUSP;}
 else
  {ib=2;}
 if((im2+NSUSP)<(n-1))
  {ie=im2+NSUSP;}
 else
  {ie=n-1;}
 for(i=ib;i<=ie;i++)
  {y[i]=yinter((ib-1),(ie+1),y[(ib-1)],y[(ie+1)],i);}
 return;
}

double dmean(void)
{
 double xm;

 if(nx==0)
 {xm=0.0;}
 else
 {xm=xs/nx;}
 return xm;
}

double stdev(void)
{
 double xsd;

 if(nx==0)
 {xsd=0.0;}
 else
 {
  if(nx==1)
  {xsd=0.0;}
  else
  {xsd=sqrt((xss-xs*xs/nx)/nx);}
 }                                /* line 352 */
 return xsd;
}

void cumul(int lmax,double x)
{
 static int ip;
 static double xq[MAXL];

 if(nx==0)
 {xs=x;
  xss=x*x;
  nx=1;
  xq[1]=x;
  ip=1;
 }
 else if(nx<lmax)
 {xs=xs+x;
  xss=xss+x*x;
  nx++;
  xq[nx]=x;
  ip=nx;
 }
 else
 {if(ip==lmax)
 {ip=1;}
 else
 {ip++;}
 xs=xs+x-xq[ip];
 xss=xss+x*x-xq[ip]*xq[ip];
 xq[ip]=x;
 }
}

double amean(int ns,int nf) 
{
 double sum_1,mean;
 int i;
/* Find the sum of matrix elements for the prescribed length */
 sum_1=0.0;
 for(i=ns;i<=nf;i++)
 {sum_1=sum_1+y[i];}
 mean=sum_1/(nf-ns);
 return mean;
}

double yinter(double x1,double x2,double y1,double y2,double xi)
{
 double yintval;
 yintval=((x2-xi)*y1+(xi-x1)*y2)/(x2-x1);
 return yintval;      /* line 402 */
}

/* Function to set up the plotting parameters */

void plotcntrl(char head[])    
{
 FILE *fout2;
 int b;
 char ch;
/* Set up plot control file used by gnuplot */
 fout2=fopen("sat_plot.gnu","w");
 b=1;
 fprintf(fout2,"set nokey\n");
 fprintf(fout2,"set data style lines\n");
 fprintf(fout2,"set title \"");
 while((b<38)&&(ch=head[b])!='_')
 {putc(ch,fout2);
  b++;}
 fprintf(fout2,"\" 0,0\n");
 fprintf(fout2,"set ylabel \"Data Value\" 0,0\n");
 fprintf(fout2,"set xlabel \"Time (sec)\" 0,0\n");
/* fprintf(fout2,"set term postscript\n"); */ /*Produce postscript file */
 fprintf(fout2,"set term pbm medium\n");
/* fprintf(fout2,"set output \"sat_plot.ps\"\n"); */
 fprintf(fout2,"set output \"output.pbm\"\n");
 fprintf(fout2,"plot \"sat_plot.dat\"\n");
 fclose(fout2);
}

/* Function to do highpass filtering */
/* Moving average filter */

void filtdath(int k)     
{
 int i;
 
 for (i=1;i<=k;i++)
 {w[i]=y[i];}
 filtdatl(k);
 for (i=1;i<=k;i++)
 {y[i]=y[i]-w[i];}
 return;
}

void filtdatl(int k)
{ 
 int m,iq,s;
 double coeff,wm,wh,c[MAXM];
 
 m=0.5*(70.0/0.1-1.0);
 coeff=1.0/(2*m+1);
 for (s=1;s<=(m+1);s++)
 {c[s]=w[s];}
 iq=1;
 wm=amean(1,(m+1));
 w[1]=wm;
 for (s=2;s<=(m+1);s++)
 {w[s]=((w[(s+m)]-wm)*coeff)+w[(s-1)];}
 for (s=(m+2);s<=(k-m);s++)
 {wh=w[s];
  w[s]=coeff*(w[s+m]-c[iq])+w[s-1];
  c[iq]=wh;
  iq=iq%(m+1)+1;}
 wm=amean((k-m),k);
 for (s=(k-m+1);s<=k;s++);
 {w[s]=coeff*(wm-c[iq])+w[s-1];
  iq=iq%(m+1)+1;}
 return;
}