External module function estates - used by mcdisp, mcdiff and singleion

This routine is optional, i.e. may be programed but is not absolutely necessary. It provides the designer of a single ion module with the possiblity to store information, usually the eigenvalues and eigenstates of a problem. This can be very useful to accelerate computations. For example, the module mcdisp provides to the single ion module functions du1calc and dMQ1 and similar at every call the matrix estates which has been set initially for every atom by a call to estates.

At the start of the programs mcdisp singleion the function estates is called (if present in the module) and it should initialize the Complex Matrix by a command such as

(*ests)=ComplexMatrix(0,nofrows,1,nofcols);

and fill this matrix with sensible numerical values for the effective field, temperature given. Parameters MODPAR and Lande Factor gJ may be used for this purpose.

The routine should look similar to

#include "vector.h"          // MatPack vector class must be included

#ifdef __linux__
extern "C" void estates(
#else
extern "C" __declspec(dllexport) void estates(
#endif
// on output
                     ComplexMatrix * ests,    // eigenstate matrix      
                                              // it must also contain population 
                                              // numbers p_i of the states (imaginary part of row 0)
                                              // and eigenvalues (real part of row 0)
// on input
                      Hxc ,                 //   vector of exchange field [meV] (can be n-dimensional, for a set of n operators)
                      Hext  ,                     // external magnetic field [T]
                      double *g_J,        // Input Lande g-factor
                      double &T,          // Input temperature (K)
                      Vector &MODPAR,     // Input vector of parameters from 
					                      //single ion property file
                      char **sipffilename)// Single ion properties filename                      
{ // ... some code to compute eigenvectors and eigenvalues

 dj=21;  // dimension matrix according to the  dimension of your eigenvalue problem
   (*ests) = ComplexMatrix(0,dj,1,dj);
   Vector En(1,dj),pn(1,dj);Matrix zr(1,dj,1,dj);Matrix zi(1,dj,1,dj);
   setup_and_solve_Hamiltonian(Hxc,Hext,En,zr,zi,sort); // call some function to solve your Hamiltonian
   
   for(i=1;i<=dj;++i)for(j=1;j<=dj;++j)(*ests)(i,j)=complex <double> (zr(i,j),zi(i,j));   // fill matrix with eigenstates
 
    //calculate partition sum and Boltzmann probabilities pn
    double zz;calculate_Z_pn(En,T,zz,pn);

     // put Energies and Boltzmann population into row 0 of eigenstates...
     for(j=1;j<=dj;++j){(*eigenstates)(0,j)=complex<double>(En(j),pn(j));}

}

Note: if $T<0$ all quantities should be evaluated assuming that all Boltzmann probabilities $p_i$ are zero except for the state number $n=(-T)$, for which the probability $p_n=1$.