/*************************************************************************/
/* Copyright (C) 1995-1997 Marlis Hochbruck, Christian Lubich, Hubert    */
/* Selhofer. All rights reserved.                                        */
/*                                                                       */
/* C code written by Mathias Froehlich.                                  */
/*                                                                       */
/* This code is part of a copyrighted package. For details, see the file */
/* "COPYING" in the top-level directory.                                 */
/*************************************************************************/

/***********************************************************************/
/* This is a time-dependent Schroedinger equation.                     */
/* It models an atom/molecule interacting with a high intensity        */
/* laser. (Peskin, Kosloff, Moiseyev: "Polynomial propagators for the  */
/* (t,t') method", J. Chem. Phys., Vol. 100, No. 12, p8849-8855.)      */
/***********************************************************************/
#include <math.h>
#include "ftypes.h"
#include "fblas.h"
#include "fft.h"
#include "exp4.h"

#define N       256
#define dx      (2.0*(a)/(N))
#define a       10.0
#define K       10.0
#define Omega   1.0
#define MU      100.0
#define T0      0.0
#define TEND    1.0
#define ATOL    1e-6
#define RTOL    1e-6

#ifndef PI
#define PI            3.14159265358979323846
#endif

extern void F77CALL (zdrscl) (finteger*,fdouble*,fdoublecomplex*,finteger*);

finteger inc1 = 1;
finteger neq = N+1;
finteger fftn = N;
fdouble fftscal = N;
fdouble constarray[N],x[N];
fdouble fftwork[4*N+15];
fdouble sqrt_K;

/* the 2 - norm of the start value, the energy of the wave package */
fdouble refenergy;

FILE *outfile;

/* time for dense output */
static fdouble dtime = T0;

void solout(long nr, fdouble told, fdouble t,fdoublecomplex *y,unsigned n,
      int *irtrn)
{
  int i;
  /* the 2 - norm of the solution, the energy of the wave package */
  fdouble energy = F77CALL (dznrm2) (&fftn,y,&inc1);

  /* print information about the actual values */
  printf("\n***************************************************************"
	"*****************\n*Solout Routine: Step number %ld\n"
	"*Time of the Problem = %e, Rel error of time variable = %e\n"
	"*Energy              = %e, Rel error of the energy    = %e\n******"
	"********************************************************************"
	"******\n\n",nr,y[fftn].r,fabs(y[fftn].r-t)/fabs(t),
	energy,fabs(energy-refenergy)/refenergy);

  /* do dense output every 1e-2 seconds */
  while (dtime < t)
  {
    fdoublecomplex work[N+1];
    ZExp4GetDenseSolution(work,1,dtime);

    fprintf(outfile,"%.16e ",dtime);
    for (i=0;i<fftn;i++)
      fprintf(outfile,"%.16e ",work[i].r);
    fputc('\n',outfile);
    dtime += 1e-2;
  }

  /* set the 'info' value to 'ok' */
  *irtrn = 0;
}

void LaserLinear(finteger sign,fdoublecomplex *fy,fdoublecomplex *y,fdouble s2)
{
  finteger i;

  s2 *= 2.0;

  F77CALL (zcopy) (&fftn,y,&inc1,fy,&inc1);

  /* forward fft */
  F77CALL (cfftf) (&fftn,fy,fftwork);

  /* do the work on the coefitients */
  for (i=0;i<fftn;i++)
  {
    fy[i].r *= constarray[i];
    fy[i].i *= constarray[i];
  }

  /* backward fft */
  F77CALL (cfftb) (&fftn,fy,fftwork);
  /* scale back ... */
  F77CALL (zdrscl) (&fftn,&fftscal,fy,&inc1);

  if (sign == 1)
    for (i=0;i<fftn;i++)
    {
      fdouble help = fy[i].r;
      fdouble fac = (K*x[i]+s2)*x[i];
      fy[i].r =   ( - fy[i].i + fac*y[i].i)/2.0;
      fy[i].i = - ( - help    + fac*y[i].r)/2.0;
    }
  else
    for (i=0;i<fftn;i++)
    {
      fdouble help = fy[i].r;
      fdouble fac = (K*x[i]+s2)*x[i];
      fy[i].r = - ( - fy[i].i + fac*y[i].i)/2.0;
      fy[i].i =   ( - help    + fac*y[i].r)/2.0;
    }    
}

void laserex(finteger neq,fdoublecomplex *fy,fdoublecomplex *y)
{
  fdouble sin_t = sin(Omega*y[neq-1].r);
  fdouble s2 = MU*sin_t*sin_t;

  LaserLinear(1,fy,y,s2);

  fy[N].r = 1.0;
  fy[N].i = 0.0;
}

void evjac(finteger neq,fdoublecomplex *res,fdoublecomplex *y,
      fdoublecomplex *vect)
{
  finteger i;
  fdouble sin_t  = sin(y[neq-1].r*Omega);
  fdouble sin_2t = sin(2.0*y[neq-1].r*Omega);
  fdoublecomplex fac;
  fdouble s2 = MU*sin_t*sin_t;

  /* the first part of the evaluation is equal to the right-hand side call */
  LaserLinear(1,res,vect,s2);


  fac.r = - vect[neq-1].i*Omega*MU*sin_2t;
  fac.i =   vect[neq-1].r*Omega*MU*sin_2t;

  for(i=0;i<neq-1;i++)
  {
    res[i].r -= x[i]*(fac.r*y[i].r - fac.i*y[i].i);
    res[i].i -= x[i]*(fac.i*y[i].r + fac.r*y[i].i);
  }

  res[neq-1].r = res[neq-1].i = 0.0;
}

int main() 
{
  finteger i;
  fdoublecomplex y0[N+1];
  fdouble t=T0;

  printf("\n#################################################################"
	"\nExp4 Example: The Laser   example, %d - dimensional Problem\n"
	"#################################################################\n"
	"\nThis example contains the use of the solout routine!\n"
	"\nPress return please ...\n",(int)neq);
  getchar();

  /* initialize fft - routines from fftpack */
  F77CALL (cffti) (&fftn,fftwork);

  /* set some constants for the rhs functions */
  sqrt_K = sqrt(K);
  for (i=0;i<(N/2);i++)
  {
    constarray[i] = ((fdouble)i)*PI/a;
    constarray[i] *= -constarray[i];
    constarray[i+N/2] = ((fdouble)(N/2-i))*PI/a;
    constarray[i+N/2] *= -constarray[i+N/2];
    x[i]     = -a + ((fdouble)i)*dx;
    x[N-i-1] = a - ((fdouble)i+1)*dx;
  }

  /* set the initial data */
  for (i=0;i<N;i++)
    y0[i].r = exp(-sqrt_K*x[i]*x[i]/2.0),y0[i].i = 0.0;

  /* the N+1-th coordinate represents the time */
  y0[N].r=t,y0[N].i=0.0;

  refenergy = F77CALL (dznrm2) (&fftn,y0,&inc1);

  /* Documentation of the functions below in ../exp4/exp4.h */

  outfile = fopen("laser.out","w");
  if (outfile == NULL) exit();

  /* define the kind of output which is used */
  ZExp4SetOutput(Solout,solout);

  /* set Atol */
  ZExp4SetAtol(Scalar,ATOL);

  /* set Rtol */
  ZExp4SetRtol(Scalar,RTOL);

  /* call the core integrator */
  ZExp4(&t,TEND,y0,neq,&laserex,NULL,&evjac,Exp4False);

  fclose(outfile);

  return 1;
}