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

int umax,ulen,ubeg,uend;
float *spc,*err,*tim,*spcx,*errx,*timx;

char *comment,*uname;


help()
{
printf("modify spectrum\n");
printf("mspc sp{[i:j]} n1 n2 ... {-sin} {-cos} {-exp} {-ln} {-sqr} {-sqrt} {-abs}\n");
printf("                         {-s n} {-d n} {-p n} {-q n} \n");
printf("the contents of sp is modified starting from i, ending with j. \n");
printf("contrary to the normal behaviour, the WHOLE spectrum (and corresponding\n");
printf("error array) is stored back \n");
printf("If n1 n2 ... is specified, then sp[i:j] is set to the constant values \n");
printf("with the last one repeated up to j. \n");
printf("other operations: \n\n");

printf("   -f   use as Filter (do not write back) \n\n");

printf("   -c text    change comment only. \n");
printf("   -tica n.m  change time calibration for this file \n");
printf("   -t    modify timespectrum only. \n");
printf("   -e    modify errorspectrum only. \n\n");

printf("   -sin  for all n: spc(n) = sin(spc(n)) \n");
printf("   -cos  for all n: spc(n) = cos(spc(n)) \n");
printf("   -exp  for all n: spc(n) = exp(spc(n)) \n");
printf("   -ln   for all n: spc(n) = log(spc(n)) \n");
printf("   -sqr  for all n: spc(n) = spc(n)*spc(n) \n");
printf("   -sqrt for all n: spc(n) = sqrt(spc(n)) \n");
printf("   -abs  for all n: spc(n) = fabs(spc(n)) \n");
printf("   -int  for all n: spc(n) = floor(spc(n)) \n");
printf("   -lin  for all n: spc(n) = n \n");
printf("   -mir  for all n: spc(n) = spc(max-n), modifies error array too \n\n");

printf("   -add x  for all n: spc(n) = spc(n) + x\n");
printf("   -sub x  for all n: spc(n) = spc(n) - x\n");
printf("   -mul x  for all n: spc(n) = spc(n) * x\n");
printf("   -div x  for all n: spc(n) = spc(n) / x\n\n");

printf("   -sr  x for all n: spc(n+x) = spc(n)  (shift right) \n");
printf("   -srt x for all n: spc(n+x) = spc(n)  (shift right, including time) \n");
printf("   -sl  x for all n: spc(n) = spc(n+x)  (shift left) \n");
printf("   -slt x for all n: spc(n) = spc(n+x)  (shift left, including time) \n");
printf("   -com x compress spectra, modifying error and timearrays too\n");
printf("   -bel x smoothly cut spectrum after x, modifies error array too\n\n");

printf("instead of numbers, other spectra arrays may be specified. \n");
printf("these spectra arrays are taken starting from the first element.\n");
printf("ONLY ONE OPERATION CAN BE PERFORMED AT A TIME ! \n");
return(0);
}

main(argc,argv)
int argc;
char *argv[];
{
int   m,n,i;
int   max;
char  *s,*z;
FILE  *stream,*fopen();


   if(argc<3) _help0(); /* at least two parameters should be specified */

   z = (char *) malloc(80);
   s = (char *) malloc(80);
   comment = (char *) malloc(80);
   uname = (char *) malloc(80);

   checkopt(argc,argv,"-dummy",z);

   spc= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   spcx= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   errx= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   timx= (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   if(timx==NULL) {
      printf("sorry, not enough memory\n");
      exit(-1);
   }

   if(checkopt(argc,argv,"-lin",z)) { /* this should be possible, without
                                         having a spectra array yet */
      for(n=1;n<=argc;n++) {
         i=instr("[",argv[n]);
         if(i>0) {
            midstr(z,argv[n],0,i-1);
         }
      }
      strcpy(s,z);
      stream=fopen(z,"r");
      if(stream==NULL) {
         strcat(z,".spc");
         stream=fopen(z,"r");
         if(stream==NULL) {            /* create a minimum spectrum */
            strcpy(z,"himmelarschundzwirn");
            writespec(s,spc,err,2,0,z);
         }
      }
   }

   
   max=readspec(argv[1],spc,err,tim,comment);
   umax=_max; ubeg=_sbeg; uend=_send;
   strcpy(uname,_spc_0nam);
   if(checkopt(argc,argv,"-o",z)) strcpy(uname,z); /* use another output */
   for(n=0;n<=umax;n++) {
      spc[n]=_uspc[n];
      err[n]=_uerr[n];
      tim[n]=_utim[n];
   }

   if(checkopt(argc,argv,"-f",z)) strcpy(uname,""); /* use as filter */
   if(checkopt(argc,argv,"-tica",z)) { /* modify time cal. */
      _tica = atosf(z);
      for(n=0;n<=umax;n++) tim[n]=_tica * n;
   }
   if(checkopt(argc,argv,"-t",z)) { /* modify time spectrum only */
      strcat(uname,".tim");
      for(n=0;n<=umax;n++) spc[n]=tim[n];
   }
   if(checkopt(argc,argv,"-e",z)) { /* modify error spectrum only */
      strcat(uname,".err");
      for(n=0;n<=umax;n++) spc[n]=err[n];
   }

   doit(argc,argv);

   writespec(uname,spc,err,umax,2,comment);
   if(instr(".tim",uname) < 0) {
      strcat(uname,".tim");
      writespec(uname,tim,err,umax,2,comment);
   }

   free(spc); free(err); free(tim); 
   free(spcx); free(errx); free(timx);
   free(uname); free(comment); free(z); free(s);
   exit(0);
}

doit(argc,argv)
int argc;
char *argv[];
{
float x,y;
int   m,n,i,j,k,max;
char  c,z[80];

   for(n=1;n<=argc;n++) {               /* check for comment modification */
      if(strcmp(argv[n],"-c")==0) {    /* modify comment */
         strcpy(comment,"");
         for(i=n+1;i<=argc;i++) {
            strcat(comment,argv[i]);
            strcat(comment," ");
         }
	 return(0);
      }
   }

   max=0;                              /* search for number or spectra name */
   for(n=2;n<argc;n++) {
      c=argv[n][0];
      if(c=='-') {
         c=argv[n][1];
         if(c>'9') break;              /* options reached, no number or file name */
         if(strlen(argv[n])==2) break;
         spcx[max++]=atof(argv[n]);     /* assume negative number */
      }
      if((c>='0') && (c<='9')) {       /* catch number(s) */
         spcx[max++]=atosf(argv[n]);     /* assume number */
	 continue;
      }
      max=readthis(argv[n],spcx,errx,timx); /* possible spectra array */
   }

   if(max>0) {                         /* set to constant */
      for(i=0;i<max;i++) spc[ubeg++]=spcx[i];
      while(ubeg <= uend) spc[ubeg++]=spcx[i-1];
      return(0);
   }

/* ----------------- binary functions ------------------ */
      
   if(checkopt(argc,argv,"-add",z)) {    /* sum */
      max=getnums(z,spcx,errx,timx);
      for(i=0;i<max;i++) {spc[ubeg] = spc[ubeg] + spcx[i]; ubeg++ ;}
      while(ubeg <= uend) {spc[ubeg] = spc[ubeg] + spcx[i-1]; ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-sub",z)) {    /* difference */
      max=getnums(z,spcx,errx,timx);
      for(i=0;i<max;i++) {spc[ubeg] = spc[ubeg] - spcx[i]; ubeg++ ;}
      while(ubeg <= uend) {spc[ubeg] = spc[ubeg] - spcx[i-1]; ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-mul",z)) {    /* product */
      max=getnums(z,spcx,errx,timx);
      for(i=0;i<max;i++) {spc[ubeg] = spc[ubeg] * spcx[i]; ubeg++ ;}
      while(ubeg <= uend) {spc[ubeg] = spc[ubeg] * spcx[i-1]; ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-div",z)) {    /* quotient */
      max=getnums(z,spcx,errx,timx);
      for(i=0;i<max;i++) {
         if(spcx[i]!=0.0) spc[ubeg] = spc[ubeg] / spcx[i];
         ubeg++ ;
      }
      while(ubeg <= uend) {
         if(spcx[i-1]!=0.0) spc[ubeg] = spc[ubeg] / spcx[i-1];
         ubeg++ ;
      }
      return(0);
   }

/* -------------------- shifting and compression ---------------- */

   if(checkopt(argc,argv,"-sr",z)) {    /* shift right */
      m=atoi(z);
      for(i=uend-m;i>=ubeg;i--) {
         spc[i+m] = spc[i];
         err[i+m] = err[i];
      }
      for(i=ubeg;i<m;i++) {
         spc[i] = 0.0;
         err[i] = 0.0;
      }
      return(0);
   }

   if(checkopt(argc,argv,"-srt",z)) {    /* shift right with time*/
      m=atoi(z);
      for(i=uend-m;i>=uend;i--) {
         spc[i+m] = spc[i];
         err[i+m] = err[i];
         tim[i+m] = tim[i];
      }
      for(i=ubeg;i<m;i++) {
         spc[i] = 0.0;
         err[i] = 0.0;
         tim[i] = 0.0;
      }
      return(0);
   }
   if(checkopt(argc,argv,"-sl",z)) {    /* shift left */
      m=atoi(z);
      for(i=ubeg;i<=uend;i++) {
         spc[i] = spc[i+m];
         err[i] = err[i+m];
      }
      for(i=uend-m;i<=uend;i++) {
         spc[i] = 0.0;
         err[i] = 0.0;
      }
      return(0);
   }
   if(checkopt(argc,argv,"-slt",z)) {    /* shift left with time */
      m=atoi(z);
      for(i=ubeg;i<=uend;i++) {
         spc[i] = spc[i+m];
         err[i] = err[i+m];
         tim[i] = tim[i+m];
      }
      for(i=uend;i<=uend;i++) {
         spc[i] = 0.0;
         err[i] = 0.0;
         tim[i] = 0.0;
      }
      return(0);
   }
   if(checkopt(argc,argv,"-com",z)) {    /* compress spectra */
      m=atoi(z);
      i=ubeg; j=ubeg;
      while(i <= (umax-m)) {
         x=0.0;
         for(k=0;k<m;k++) x=x+spc[i++];
         spc[j] = x/m;
         err[j] = err[i] / sqrt((double)m);
         tim[j] = tim[j] * m;
         j = j + 1;
      }
      umax=j-1;
      return(0);
   }

   if(checkopt(argc,argv,"-mir",z)) {    /* mirror spectra */
      m=atoi(z);
      i=ubeg; j=umax;
      while(i < j) {
         x=spc[i];
         spc[i]=spc[j];
         spc[j]=x;
         x=err[i];
         err[i]=err[j];
         err[j]=x;
         j=j-1; i=i+1;
      }
      return(0);
   }

   if(checkopt(argc,argv,"-bel",z)) {    /* cut spectrum smoothly*/
      m=atoi(z);
      i=ubeg; j=uend;
      while(i < j) {
         x =  ((float) i) / ((float) m);
         x=cos(x*1.57079);
         if(i>=m) x=0.0;
         spc[i] = x * spc[i];
         err[i] = x * err[i];
         i=i+1;
      }
      return(0);
   }
         
/* ----------------------- monaric functions ------------------ */

   if(checkopt(argc,argv,"-sin",z)) {
      while(ubeg <= uend) {spc[ubeg] = sin(spc[ubeg]); ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-cos",z)) {
      while(ubeg <= uend) {spc[ubeg] = cos(spc[ubeg]); ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-exp",z)) {
      while(ubeg <= uend) {spc[ubeg] = exp(spc[ubeg]); ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-ln",z)) {
      while(ubeg <= uend) {spc[ubeg] = log(spc[ubeg]); ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-sqr",z)) {
      while(ubeg <= uend) {spc[ubeg] = spc[ubeg] * spc[ubeg]; ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-sqrt",z)) {
      while(ubeg <= uend) {spc[ubeg] = sqrt(spc[ubeg]); ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-abs",z)) {
      while(ubeg <= uend) {spc[ubeg] = fabs(spc[ubeg]); ubeg++ ;}
      return(0);
   }
   if(checkopt(argc,argv,"-int",z)) {
      while(ubeg <= uend) {spc[ubeg] = floor(spc[ubeg]); ubeg++ ;}
      return(0);
   }
      
}

readthis(s,spc,err,tim)
char s[];
float spc[],err[],tim[];
{
int   max;

   max=readspec(s,spc,err,tim,s);
   return(max);
}

getnums(z,spc,err,tim)
char z[];
float spc[],err[],tim[];
{
int i,n,max;

   max=0;
   if((z[0]=='-') || (z[0]<='9')) {
      spc[max++]=atof(z);
   } else {
      max=readthis(z,spc,err,tim); /* possible spectra array */
   }
   return(max);
}