External module function du1calc - used by mcdisp, singleion

The external single ion module has to provide the components nofcomponentsxnofcomponents matrix $M^s_{\alpha\beta}$ (see equation (230)) for every transition $\vert-\rangle \rightarrow \vert+\rangle$which is to be taken into consideration in the calculation. Note, in general $M^s_{\alpha\beta}$ it is a quadratic matrix with the same dimension as the vectors $\mathbf I$ and Hxc. Actually it is not necessary to compute the whole Matrix $M^s_{\alpha\beta}$, because only the first eigenvector is needed in the computations. Thus, in order to make calculations easier and provide a unique phase of the eigenvectors, the external single ion module must return not the Matrix $M$ but the unnormalized eigenvector $u1$, which is given by equation (232): ${\mathbf u^s_{\alpha1}}=\sqrt{(p_--p_+)}\langle -\vert I^s_{\alpha}-\langle I^s_{\alpha}\rangle_{\mathbf H,T}\vert+\rangle$. Note that in contrast to ${\mathcal U^s_{\alpha1}}$ the eigenvector ${\mathbf u^s_{\alpha1}}$ is not normalised (no $\gamma$!) and thus the matrix $M^s_{\alpha\beta}$ may be recovered from it.

The format to be used is:

extern "C" int du1calc(int & tn,double & T,Vector & Hxc,Vector & Hext,double * g_J,Vector & MODPAR,
char ** sipffilename,ComplexVector & u1,float & delta, int & n, int & nd, ComplexMatrix & est)

The meaning of the symbols is as follows:

on input
   |tn|            transition-number  
   sign(tn)        >0 standard, <0 routine should do some printout to stdout for user information
   MODPAR          Vector with Parameters  read in single ion property file
   sipffilename    file name of the single ion parameter file
   g_J             Lande factor
   T               Temperature[K]
   Hxc             vector of exchange field [meV] (can be n-dimensional, for a set of n operators)
   Hext            external magnetic field [T]
   est             eigenstate matrix (initialized by estates)
                   it should/may also contain population numbers of the states
				   (imaginary part of row 0)
                   and eigenvalues (real part of row 0) with values set by the most recent call
				   for this ion (use of this matrix is optional)
   u1(1)           ninit + i pinit (from mcdisp options  -ninit and -pinit)
   delta           upper boundary for transition energy (meV) to be considered
                   (from mcdisp option -maxE)
on output
   int             total number of transitions
   delta           transition energy [meV]
   n,nd            initial and final level
   u1             vector u1=<-|Ialpha-<Ialpha>|+>sqrt((p- - p+))                          
                note that as in Icalc the single ion momentum vector <-|I|+> 
				(n- dimensional with n>=1)
                may be an arbitrary set of operators

The module function must perform the following tasks:

1. check if the dimensions of vectors I,Hxc (taken by mcphas from the number of interaction constant columns in mcphas.j) and MODPAR (taken by mcphas from the number of params in the single ion property file) agree with the module specifications. Note a module may be designed to take different dimensions depending on the input files, however the dimensions of vectors I and Hxc has to agree and must be within the range of dimensions which can be treated by the module. If the check fails the module function should exit the program with an appropriate error message
2. the module function should do a numbering of all possible single ion transitions and return the total number of transitions as an integer. Input file parameters params are supplied as a vector MODPAR and Lande factor as g_J and can be used for this purpose. The numbering will depend on the parameters ninit, pinit and maxE which are provided as input. These parameters have to be considered. IMPORTANT: the numbering scheme of transitions has to be the same for du1calc and all the corresponding d...1 functions for observables !
3. it should calculate from meanfields at a given temperature the transition energy of transition number tn. The result should be returned as delta
4. for the transition number tn the vector u1 is to be filled with ${\mathbf u^s_{\alpha1}}=\sqrt{(p_--p_+)}\langle -\vert\hat I^s_{\alpha}-\langle \hat I^s_{\alpha}\rangle_{\mathbf H,T}\vert+\rangle$.
5. If the energy of this transition is zero, i.e. $\Delta(tn)=0$ (diffuse scattering), the expression (230) would be zero because $(p_--p_+)$ vanishes. In this case the single ion module should calculate $(p_+/kT)$ instead of $(p_--p_+)$.
6. if transition is from a level to itself, then a negative value of $\Delta=-10^{-10}$ should be returned.
7. 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$.