/*
	test harness for optimisers: 
	uses variational principle to find ground state of Schroedinger's Hamiltonian
*/

#include<iostream>
using namespace std;

#include<wxmaths.h>
using namespace wxmaths;

class Complex_vector: public Array<Complex> { };

class wavefunction: public Optimisable_function
{
   public:
   
   double hbar, m, q; // Planck's constant / 2 pi; quantum mass, charge
 
   Array<double> vol; // volume
   Array<double> v; // electrostatic potential
   Array<Complex_vector> phi; // basis functions
   Array<Complex_vector> d2phi; // basis functions' Laplacians
   
   Array<double> amp; // basis function amplitudes
   
   void Read(istream& s)
   {
       s >> hbar >> m >> q;
      int npts, nfns; s >> npts >> nfns;
      vol.Recreate(npts); v.Recreate(npts);
      phi.Recreate(npts); d2phi.Recreate(npts);
      vol.Read_elements(s); v.Read_elements(s);
      for (int i=0;i<npts;i++) {
         phi[i].Recreate(nfns); d2phi[i].Recreate(nfns);
         phi[i].Read_elements(s); d2phi[i].Read_elements(s);
      }
      
      amp.Recreate(nfns); amp.Read_elements(s);
   }
   void Write(ostream& s) const
   {
      s << hbar << ' ' << m << ' ' << q << endl;
      int npts = v.n, nfns = phi[0].n; s << npts << ' ' << nfns << endl;
      vol.Write_elements(s); s << endl;
      v.Write_elements(s); s << endl;
      for (int i=0;i<npts;i++) {
         phi[i].Write_elements(s); d2phi[i].Write_elements(s); s << endl;
      }
      amp.Write_elements(s); s << endl;
   }
   
   double number(Vector& c) const
   {
      double cum = 0.0;
      for (int i=0;i<vol.n;i++) {
         Complex psi(0.0,0.0);
         for (int j=0;j<phi[i].n;j++) psi = psi + c[j] * (phi[i])[j];
         cum = cum + vol[i] * norm(psi);
      }
      return cum;
   }
   
   int constrain(Vector& c) const
   {
      double wvint = ::sqrt(number(c));
      double fac = 1.0 / wvint;
      for (int i=0;i<c.n;i++) c[i] = c[i] * fac;
      return 1;
   }
   
   double evaluate(const Vector& c) const
   {
      double h = 0.0;
      for (int i=0;i<vol.n;i++) {
         Complex psi(0.0,0.0), d2psi(0.0,0.0);
         for (int j=0;j<phi[i].n;j++) {
            psi = psi + c[j] * (phi[i])[j];
            d2psi =d2psi + c[j] * (d2phi[i])[j];
         }
         h = h + real(vol[i] * conjugate(psi) * (v[i] * psi - (hbar / (2 * m)) * d2psi));
      }
      return h;
   }
};

int main()
{
wavefunction w; cin >> w;
Optimiser *o; cin >> o;
Vector c(w.amp); o->minimise(w,c);
cout << "final values: " << c << endl;

return 0;
}