/*
   ratio function:
   sp = (usp1 - usp2) / (usp1 + usp2)
*/

#include <stdio.h>
#include <spec.h>
#define TMPFILE "global.def"

float *spc, *err, *tim, *ZNu, *ZSu, *ZNd, *ZSd, *A, *au, *ad;
int   flg_v;

help()
{
  printf("\n\nRatio Function:  au = (ZNu - ZSu) / (ZNu + ZSu)  ; u = up\n");
  printf("                 ad = (ZNd - ZSd) / (ZSd + ZSd)  ; d = down\n");
  printf("                 A = 0.5 * (ad - au)\n");
  printf("Error:           dax = 1/sqrt(ZNx + ZSx)  ;  x = u,d\n");
  printf("                 A = 0.5 * sqrt((dad * dad) + (dau * dau))\n");
  printf("nmrrvt filename [-ou file] [-od file] [-oA file]\n");
  printf("       filename must be specified without numbers and extension.\n");
  printf("       filename1.spc is expected to be ZNu (North Up)\n");
  printf("       filename2.spc is expected to be ZSu (South Up)\n");
  printf("       filename3.spc is expected to be ZNd (North Down)\n");
  printf("       filename4.spc is expected to be ZSd (South Down)\n");
  printf("options:\n");
  printf("   -ou filename  output file for UP instead of /tmp/au.spc\n");
  printf("   -od filename  output file for DOWN instead of /tmp/ad.spc\n");
  printf("   -oA filename  output file for 0.5(ad - au) instead of /tmp/A.spc\n");
  printf("\n(C) Rainer Kowallik\n\n");
  exit(0);
}

float sd(x,y)
float x,y;
{
if(y!=0.0) return(x/y);
return(0.0);
}

main(argc,argv)
int argc;
char *argv[];
{
int  n,m,i,max,flg_down,lastn,xbeg,xend;
char z[80],comment[80],nam_u[80],nam_d[80],nam_a[80];
FILE *fp;

   strcpy(nam_u,"/tmp/au");
   strcpy(nam_d,"/tmp/ad");
   strcpy(nam_a,"/tmp/A");
   flg_down = 0;            /* defaults to "down spectra absent" */

   if(checkopt(argc,argv,"-ou",z)) strcpy(nam_u,z); 
   if(checkopt(argc,argv,"-od",z)) strcpy(nam_d,z); 
   if(checkopt(argc,argv,"-oA",z)) strcpy(nam_a,z);

/* -----------------------------------------------
   reading spectrum 
   ----------------------------------------------- */

   spc=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   err=(float *)calloc(_MAXSPCLEN+2,sizeof(float));
   tim=(float *)calloc(_MAXSPCLEN+2,sizeof(float));

   sprintf(z,"%s1",argv[1]);
   max=readspec(z,spc,err,tim,comment); /* now we know the spectra length ! */

   ZNu = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   ZSu = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   ZNd = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   ZSd = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   A   = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   au  = (float *)calloc(_MAXSPCLEN+2,sizeof(float));
   ad  = (float *)calloc(_MAXSPCLEN+2,sizeof(float));

   for(n = 0; n < max; n++) ZNu[n] = spc[n];
   sprintf(z,"%s2",argv[1]);
   max=readspec(z,ZSu,err,tim,comment);

   sprintf(z,"%s3.spc",argv[1]);
   fp = fopen(z,"r"); /* check existance of "down" spectra */
   if(fp != NULL) {
      flg_down = 1;
      fclose(fp);
   }
   if(flg_down != 0) {
      sprintf(z,"%s3",argv[1]);
      max=readspec(z,ZNd,err,tim,comment);
      sprintf(z,"%s4",argv[1]);
      max=readspec(z,ZSd,err,tim,comment);
   }


/* -------------------------------------------
           calculations
   ------------------------------------------- */
   for(n = 0; n < max; n ++) {
      au[n] = sd((ZNu[n] - ZSu[n]) , (ZNu[n] + ZSu[n]));
      err[n] = sd(1.0,(float)sqrt(ZNu[n] + ZSu[n]));
   }
   writespec(nam_u,au,err,max,2,comment);
   strcat(nam_u,".tim");
   writespec(nam_u,tim,err,max,2,comment);


   if(flg_down != 0) {
      for(n = 0; n < max; n ++) {
         ad[n] = sd((ZNd[n] - ZSd[n]) , (ZNd[n] + ZSd[n]));
         err[n] = sd(1.0,(float)sqrt(ZNd[n] + ZSd[n]));
      }
      writespec(nam_d,ad,err,max,2,comment);
      strcat(nam_d,".tim");
      writespec(nam_d,tim,err,max,2,comment);
      for(n = 0; n < max; n ++) {
         A[n] = 100.0 * (ad[n] - au[n]);
         err[n] = 100.0 * 
  (float)sqrt((double)(sd(1.0,ZNu[n] + ZSu[n]) + sd(1.0,ZNd[n] + ZSd[n])));
      }
      writespec(nam_a,A,err,max,2,comment);
      strcat(nam_a,".tim");
      writespec(nam_a,tim,err,max,2,comment);
   }
}