/********************************************************************/
/*                                                                  */
/*    MODUL              : SCALC.C                                  */
/*    VERSION            : 1.0                                      */
/*    DATE               : 24/09/91                                 */
/*    WRITTEN BY         : Christoph Bohle                          */
/*    Modified for spectra calculation by Rainer Kowallik, 06/11/91 */
/* ---------------------------------------------------------------- */
/*    LANGUAGE           : C                                        */
/*    OS                 : UNIX                                     */
/*                                                                  */
/********************************************************************/

/********************************************************************/
/*   INCLUDE-Files                                                  */
/********************************************************************/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include <spec.h>

#ifndef RAND_MAX
#define RAND_MAX 2147483647
#endif
 
int spclenmax = 0;

help()
{
printf("scalc  spectra calculator. Call convention:\n");
printf("    scalc [newspc =] formula(function .. number .. spectrum ..)\n");
printf("\nIf no output spectrum is specified, only the first element is\n");
printf("printed on standard output.\n");
printf("Outputspectrum can be specified with a range.\n");
printf("If outputspec allready exists, only the specified part will\n");
printf("be substituted by the result of the calculation.\n");
printf("No calculation will be done with the error and time arrays !\n");
printf("\nPossible function:\n");
printf("exp(expression)\n");
printf("log(expression)\n");
printf("log10(expression)\n");
printf("sin(expression)\n");
printf("cos(expression)\n");
printf("tan(expression)\n");
printf("asin(expression)\n");
printf("acos(expression)\n");
printf("atan(expression)\n");
printf("sinh(expression)\n");
printf("cosh(expression)\n");
printf("tanh(expression)\n");
printf("sqrt(expression)\n");
printf("rnd()                  returns a random number between 0 and 1\n");
printf("lin()                  for(n=1;n<_MAXSPCLEN;n++) spc[n] = n;\n");
printf("sum(spc)               for(n=1;n<_MAXSPCLEN;n++) sum = sum + spc[n];\n");
printf("pi()                   for(n=1;n<_MAXSPCLEN;n++) spc[n] = 3.1415927;\n");
printf("\nOperators are +,-,*,/ and ^\n");
printf("   programmed by Christoph Bohle and Rainer Kowallik\n");
exit(0);
}


       
/********************************************************************/
/*    Functions                                                     */
/********************************************************************/


float *expression();

char *end_of_brackets( str )
char *str;  
{                                   
  int i;
  i=1;
  do
  {	
    str++;
    switch( *str )
    {
      case '('  : i++;
                  break;
      case ')'  : i--;
	  	            break;
      case '\0' : i=0;
    }
  }
  while( i>0 );
  return( str );
}
		
float *factor( facstr )
char *facstr;
{
int i,n;
float *result, *tim, *err;
float x,y;
char *help, c;
char comment[80];
  
  /* Float */
  if( isdigit(*facstr) || *facstr=='.' || *facstr=='-' || *facstr=='+' ) {
    result = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
    y = atof(facstr);
    for(i = 0; i < _MAXSPCLEN; i++) result[i] = y;
    return(result);
  }
  
  /* Expression in brackets */
  if( *facstr=='(' )
  {
    *(end_of_brackets( ++facstr )) = '\0';
    return( expression( facstr ) );
  }
  
  /* Function */
  help = facstr;
  while( isalnum(*help) || *help == '_' || *help == '.'
                        || *help == '[' || *help == ':' || *help == ']') help++;
  if( help != facstr )
  {
    c = *help;
    *help = '\0';
    
    /* Function */
    if( strcmp( facstr , "exp" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = exp(result[i]);
      return( result );
    }
    if( strcmp( facstr , "log" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = log(result[i]);
      return( result );
    }
    if( strcmp( facstr , "log10" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = log10(result[i]);
      return( result );
    }
    if( strcmp( facstr , "sin" ) == 0 )
    {
      *help = c;      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = sin(result[i]);
      return( result );
    }
    if( strcmp( facstr , "cos" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = cos(result[i]);
      return( result );
    }
    if( strcmp( facstr , "tan" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = tan(result[i]);
      return( result );
    }
    if( strcmp( facstr , "asin" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = asin(result[i]);
      return( result );
    }
    if( strcmp( facstr , "acos" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = acos(result[i]);
      return( result );
    }
    if( strcmp( facstr , "atan" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = atan(result[i]);
      return( result );
    }
    if( strcmp( facstr , "sinh" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = sinh(result[i]);
      return( result );
    }
    if( strcmp( facstr , "cosh" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = cosh(result[i]);
      return( result );
    }
    if( strcmp( facstr , "tanh" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = tanh(result[i]);
      return( result );
    }
    if( strcmp( facstr , "sqrt" ) == 0 )
    {
      *help = c;
      result = factor( help );
      for(i = 0; i < spclenmax; i++) result[i] = sqrt(result[i]);
      return( result );
    }
    if( strcmp( facstr , "rnd" ) == 0 )            /* generate random array */
    {
      *help = c;
      result = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
      for(i = 0; i < _MAXSPCLEN; i++) result[i] = ((float) rand() / (RAND_MAX + 1.0)) ;
      return( result );
    }
    if( strcmp( facstr , "lin" ) == 0 )            /* generate linear array */
    {
      *help = c;
      result = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
      for(i = 0; i < _MAXSPCLEN; i++) result[i] = (float) i ;
      return( result );
    }
    if( strcmp( facstr , "pi" ) == 0 )
    {
      *help = c;
      result = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
      for(i = 0; i < _MAXSPCLEN; i++) result[i] = 3.1415927 ;
      return( result );
    }
    if( strcmp( facstr , "sum" ) == 0 )
    {
      *help = c;
      result = factor( help );
      y = 0.0;
      for(i = 0; i < spclenmax; i++) y = y + result[i];
      for(i = 0; i < spclenmax; i++) result[i] = y;
      return( result );
    }

/*  read spectra array */
    result = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
    tim = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
    err = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
    i = readspec(facstr,result,err,tim,comment);
    if(i > spclenmax) spclenmax = i;
    free(err); free(tim);
    return( result );
  }

  /* Default */
  if( *facstr == '\0' ) {
    result = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
    for(i = 0; i < _MAXSPCLEN; i++) result[i] = 0.0;
    return( result );
  } else {
    return( factor( facstr+1 ) );
  }
}

float *po( postr )
char postr[];
{
char *help;
float *ar1, *ar2;
int i,n;

  help = postr;
  while( (*help != '^') && (*help != '\0' ) )
  {
    if( *help == '(' )
    	help = end_of_brackets( help );
    help++;
  }
  
  switch( *help )
  {
	  case '^' : *help = '\0';
               ar1 = factor( postr );
               ar2 = po( help + 1);
               for(i = 0; i < spclenmax; i++) ar1[i] = pow(ar1[i],ar2[i]);
               free(ar2);
               return( ar1 );
    case '\0' : if( help == postr ) {
                  ar1 = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
                  for(i = 0; i < spclenmax; i++) ar1[i] = 0.0;
                  return( ar1 );
                } else {
                  return( factor( postr ) );
                }
  }
}

float *term( termstr )
char termstr[];
{
char *help;
float *ar1, *ar2;
int i,n;
  
  help = termstr;
  while( (*help != '*') && (*help != '/') && (*help != '\0' ) )
  {
    if( *help == '(' )
    	help = end_of_brackets( help );
    help++;
  }
  
  switch( *help )
  {
	 case '*' : *help = '\0';
               ar1 = po( termstr );
               ar2 = term( help + 1);
               for(i = 0; i < spclenmax; i++) ar1[i] = ar1[i] * ar2[i];
               free(ar2);
               return( ar1 );
    case '/' : *help = '\0';
               ar1 = po( termstr );
               ar2 = term( help + 1);
               for(i = 0; i < spclenmax; i++) {
                  if(ar2[i] != 0.0) ar1[i] = ar1[i] / ar2[i];
               }
               free(ar2);
               return( ar1 );
    case '\0' : if( help == termstr ) {
                  ar1 = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
                  for(i = 0; i < spclenmax; i++) ar1[i] = 0.0;
                  return( ar1 );
                } else {
                  return( po( termstr ) );
                }
  }
}

float *expression( exprstr )
char exprstr[];
{
char *help;
float *ar1, *ar2;
int i,n;

  help = exprstr;
  while( (*help != '+') && (*help != '-') && (*help != '\0' ) )
  {
    if( *help == '(' )
    	help = end_of_brackets( help );
    help++;
  }
  
  switch( *help )
  {
	 case '+' : *help = '\0';
               ar1 = term( exprstr );
               ar2 = expression( help + 1);
               for(i = 0; i < spclenmax; i++) ar1[i] = ar1[i] + ar2[i];
               free(ar2);
               return( ar1 );
    case '-' : *help = '\0';
               ar1 = term( exprstr );
               ar2 = expression( help + 1);
               for(i = 0; i < spclenmax; i++) ar1[i] = ar1[i] - ar2[i];
               free(ar2);
               return( ar1 );
    case '\0' : if( help == exprstr ) {
                  ar1 = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
                  for(i = 0; i < spclenmax; i++) ar1[i] = 0.0;
                  return( ar1 );
                } else {
                  return( term( exprstr ) );
                }
  }
}

float *analyse_function( funstr )
char *funstr;
{
char *help;
float *result;
  
  help=malloc( strlen(funstr)+3 );
  strcpy( help , funstr );
  result = expression( help );
  free( help );
  return( result );  
}

/********************************************************************/
/* Main-function                                                    */
/********************************************************************/

main( argc , argv  )
int argc;
char *argv[];              
{
int i, l, ix1, ix2;
char *fun, *snam;
float *result, *oldspc, *err;
  
  snam = (char *) malloc(80);

  i = checkopt(argc,argv,"-sowatjibtsnich",snam); /* we must call checkopt once */

  l = 1;
  for( i=1 ; i < argc ; i++ )
  {
    l += strlen( argv[i] ) + 1;
  }
  fun = malloc( l + 3 ); *fun = '\0';
  for( i=1 ; i < argc ; i++ )
  {
    strcat( fun , argv[i] );
    strcat( fun , " ");
  }

  ix1 = 0;
  ix2 = 0;
  oldspc = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
  err = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
  find_res_spc(fun,snam,&ix1,&ix2,oldspc,err);
  result = analyse_function( fun );

  if((ix1 != 0) || (ix2 != 0)) {
     if(ix2 < 0) ix2 = spclenmax;
     for(i = ix1; i <= ix2; i++) oldspc[i] = result[i];
     writespec(snam,oldspc,err,spclenmax,2,"test");
     free(err); free(oldspc);
  } else {
     printf("%f\n",(double)result[0]);
  }
  free(snam); free(result);
  exit(0);
}

find_res_spc(str,snam,ix1,ix2,oldspc,err)
char *str, *snam;
int  *ix1, *ix2;
float *oldspc, *err;
{
int i,j,l,n;
char c;
char *z, *s;
FILE *stream;
float *tim;


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

   l = instr("=",str);
   if(l < 0) return;

   snam[0] = 0; i = 1; j = 0; c = str[0];
   while((c != 0) && (c != '[') && (c != '=')) {
      if(c != ' ') snam[j++] = c;
      c = str[i++];
   }
   snam[j] = 0;

   if(c == '[') {
      s[0] = 0; j = 0; c = str[i++];
      while((c != 0) && (c != ':') && (c != ']')) {
         if(c != ' ') s[j++] = c;
         c = str[i++];
      }
      s[j] = 0; *ix1 = atoi(s); *ix2 = *ix1;
      if(c == ':') {
         s[0] = 0; j = 0; c = str[i++];
         while((c != 0) && (c != ']')) {
            if(c != ' ') s[j++] = c;
            c = str[i++];
         }
         s[j] = 0; *ix2 = atoi(s);
      }
   } else {
      *ix1 = 0;
      *ix2 = -1;
   }

   j = 0; n = l; l = strlen(str);
   for(i = n; i <= l; i++) str[j++] = str[i];
  
   strcpy(z,snam);
   stream = fopen(z,"r");
   if(stream == NULL) {
      strcat(z,".spc");
      stream=fopen(z,"r");
   }
   if(stream!=NULL) {
      fclose(stream);
      tim = (float *) calloc(_MAXSPCLEN+2,sizeof(float));
      i = readspec(snam,oldspc,err,tim,s);
      free(tim);
      if(i > spclenmax) spclenmax = i;
   }
   
   if(*ix2 > spclenmax) spclenmax = *ix2;
   free(s); free(z);
   return(0);
}