// derivs.cc

#include "RKtest.h"

void derivs(const double x, const int nvar, double *V, double *dVdx)
  // This routine computes the diff. eqn. for the Burgers equation 
  // in spectral form.  It must be used only with spectra that are 
  // derived by taking:
  //             N
  //        P = SUM P_n exp(-jny)
  //            n=1
  // It assumes that P_n= Re(P_n)+j*Im(P_n) and that A=A_n-j*Delta_n 
  // If the sum has a one half in front of it, then the "Factor of n/2"
  // below should be changed to a "Factor of n/4".
{
  int nhar;                 // Number of harmonics
  int dm, dnMm, dmMn;       // Temporary variables
  double tr, ti;            // Real and imag parts of deriv (temporary)
  
  nhar=(int)(nvar/2);
  for (int n=1;n<=nhar;n++){                 // Compute for all harmonics
    tr=ti=0;                                 // Initialize variables
    for (int m=n+1;m<=nhar;m++){             // Sum from n+1 to nhar
      dm=2*m;
      dmMn=2*(m-n);
      tr+= -V[dm-1]*V[dmMn  ]+V[dm]*V[dmMn-1];
      ti+= -V[dm-1]*V[dmMn-1]-V[dm]*V[dmMn  ];
    }
    tr*=2.0;                                 // Factor of two
    ti*=2.0;
    for (int m=1;m<=n-1;m++){                // Sum from 1 to n-1
      dm=2*m;
      dnMm=2*(n-m);
      tr+=  V[dm-1]*V[dnMm  ]+V[dm]*V[dnMm-1];
      ti+= -V[dm-1]*V[dnMm-1]+V[dm]*V[dnMm  ];
    }
    tr*=n/2.0;                               // Factor of n/2
    ti*=n/2.0;
    //tr*=n/4.0;                             // Factor of n/4
    //ti*=n/4.0;
    
    tr+= -n*n*A*V[2*n-1];                    // Adding dissipation
    ti+= -n*n*A*V[2*n  ];

    tr+= -(-n*kxbar*(n*D/(1+n*D)))*V[2*n  ]; // Adding dispersion
    ti+=  (-n*kxbar*(n*D/(1+n*D)))*V[2*n-1];

    dVdx[2*n-1]=tr;                          // Final result
    dVdx[2*n  ]=ti;
    
  }
  
}