// CalcB.cc
// 
// Original author:  Yurii Il'inskii
//
// Modifications:
// 08 Dec 1996: Modified by Ronald Kumon; realigned text and added comments;
//              included introductory comments.  
// 07 Feb 1998: Added routine ShadeHarmonics to provide a shading function
//              during the waveform reconstruction from the harmonics.

#include <math.h>
#include <string.h>
#include <iostream.h>
#include "Calc.h"

const double Pi=3.141592653589793;

void ShadeHarmonics(int nm,double *V,double *Vshade)
// Shades harmonic amplitudes by a function specified within.
{
  int nm2;          // Total number of harmonic components
  double mid;       // Parameter that defines width of shading
  double exponent;  // Exponent for shading function

  mid=70.0;
	
  for (int n=1;n<=nm;n++) {
    //exponent=pow(n/mid,16);
    // cout << "n= " << n << " Factor= " << exp(-exponent) << endl;
    //Vshade[2*n-1]=V[2*n-1]*exp(-exponent);
    //Vshade[2*n  ]=V[2*n  ]*exp(-exponent);
    Vshade[2*n-1]=V[2*n-1];
    Vshade[2*n  ]=V[2*n  ];
  }
}

void CalcVx(int nm,double *V,double BxrN,double BxiN,double *Vx)
// Computes velocity waveform in x-direction 
// given harmonic amplitudes
{
  int nm2;        // Total number of harmonic components
  double dT,T;    // Step size variables
	
  dT=Pi/nm;       // Fourier step size 
  nm2=nm*2;       // Real and imag components cause doubling
	
  for (int n=0;n<=nm2;n++) { // Compute Vx[k]=Sum Re(V[n]*exp(-i*Pi*k*n/nm))
    Vx[n]=0.0;
    for (int m=1;m<=nm;m++) {
      T=dT*m*n;
      Vx[n] += V[2*m-1]*cos(T)+V[2*m]*sin(T);
    }
    Vx[n]*=2.0;   // Factor from adding term and conjugate
  }
}

void CalcVxh(int nm,double *V,double BxrN,double BxiN,double *Vxh)
// Computes Hilbert transformed velocity waveform in x-direction 
// given harmonic amplitudes
{
  int nm2;       // Total number of harmonic components
  double dT,T;   // Step size variables
	
  dT=Pi/nm;      // Fourier step size 
  nm2=nm*2;      // Real and imag components cause doubling
	
  for (int n=0;n<=nm2;n++) { // Compute Vx[k]=Sum Re(i*Bx*V[n]*exp(-i*Pi*k*n/nm))
    Vxh[n]=0.0;
    for (int m=1;m<=nm;m++) {
      T=dT*m*n;
      Vxh[n] +=  (BxrN*V[2*m-1]-BxiN*V[2*m])*sin(T)
                -(BxrN*V[2*m]+BxiN*V[2*m-1])*cos(T);
    }
    Vxh[n]*=2.0; // Factor from adding term and conjugate
  }
}

void CalcVy(int nm,double *V,double ByrN,double ByiN,double *Vy)
// Computes velocity waveform in y-direction 
// given harmonic amplitudes
{
  int nm2;       // Total number of harmonic components
  double dT,T;   // Step size variables
	
  dT=Pi/nm;      // Fourier step size 
  nm2=nm*2;      // Real and imag components cause doubling
	
  for (int n=0;n<=nm2;n++) { // Compute Vy[k]=Sum Re(By*V[n]*exp(-i*Pi*k*n/nm))
    Vy[n]=0.0;
    for (int m=1;m<=nm;m++) {
      T=dT*m*n;
      Vy[n] += (ByrN*V[2*m-1]-ByiN*V[2*m])*cos(T)
              +(ByrN*V[2*m]+ByiN*V[2*m-1])*sin(T);
    }
    Vy[n]*=2.0;  // Factor from adding term and conjugate
  }	
}

void CalcVz(int nm,double *V,double BzrN,double BziN,double *Vz)
// Computes velocity waveform in z-direction 
// given harmonic amplitudes
{
  int nm2;       // Total number of harmonic components
  double dT,T;   // Step size variables
	
  dT=Pi/nm;      // Fourier step size 
  nm2=nm*2;      // Real and imag components cause doubling
	
  for (int n=0;n<=nm2;n++) { // Compute Vz[k]=Sum Re(Bz*V[n]*exp(-i*Pi*k*n/nm))
    Vz[n]=0.0;
    for (int m=1;m<=nm;m++) {
      T=dT*m*n;
      Vz[n] += (BzrN*V[2*m-1]-BziN*V[2*m])*cos(T)
              +(BzrN*V[2*m]+BziN*V[2*m-1])*sin(T);
    }
    Vz[n]*=2.0;  // Factor from adding term and conjugate
  }
}

void createSuffix(const char *datastr, char *newstr)
// This function takes a string that contains a number, prepends a 
// prefix and postpends a suffix and returns the new string.  It is
// used to make the suffixes for the position data files.
{
  const char prefix[]="-";
  strcpy(newstr,prefix);
  strcat(newstr,datastr);
}