/* standard begin
*/

#define MAXTERM 30

#include <stdio.h>
#ifdef AMIGA
#include <gfxamiga.h>
#include <libraries/dos.h>
#include <libraries/dosextens.h>
#else
#include <gfx.h>
#endif

float  *spc,*err,*tim;

double x,y;
double Pi = 3.1415927,
       Gs = 2.0,
       k = 1.38062E-23,
       Mu0 = 1.256637E-6,
       MuB = 9.2741E-24,
       gI = 0.06,
       Nq = 4.0,
       l = 3.0,
       Nel = 5.0,
       J_grund = 0.0,
       Fak_BN = 0.0,  /* must be set within programm */
       R_mittel = 0.0, /* must be set within programm */
       B_core = 10.0,
       Tbeg = 50.0,
       Tend = 1000.0,
       Tstep = 10.0;

double *EJf,
       *Lges,
       *Sges,
       *Jges;

double *sum1,
       *sum2,
       *sum3;

double CurieTerm,
       VanVleckTerm;

int    NTerm = 0,
       Flag_1dT = 1;

char   *fnam;

help()
{
  printf("this programm calculates the paramagnetic enhancement factor for\n");
  printf("various ions at various temperatures.\n");
  printf("B_HF [-r file] [-T] [-Ts n.m] [-Te n.m] [-Ti n.m] [-go]\n");
  printf("      -r file    reads data from file\n");
  printf("      -T         generates spectrum as function of direct Temperature\n");
  printf("      -Ts n.m    sets start temperature\n");
  printf("      -Te n.m    sets end temperature\n");
  printf("      -Ti n.m    sets temperature increment\n");
  printf("      -go        suppress interactive mode\n");
}

double dabs(x)
double x;
{
   if(x >= 0.0) return(x);
   return(-1.0 * x);
}

double Energy(E,us,dir)
double E;
char *us;
int dir;
{
double c;

   c = 0.0;
   if(strcmp(us,"J") == 0) c = 1.0;
   if(strcmp(us,"K") == 0) c = 1.38062E-23;
   if(strcmp(us,"eV") == 0) c = 1.60219E-19;
   if(strcmp(us,"erg") == 0) c = 1.0E-7;
   if(c == 0.0) {
      fprintf(stderr,"unrecognised Energy unit :%s\n",us);
      return(E);
   }
   if(dir < 0) c = 1/c;
   return(c * E);
}

double ReCuR(R,us,dir)
double R;
char *us;
int dir;
{
double c;

   c = 0.0;
   if(strcmp(us,"m") == 0) c = 1.0;
   if(strcmp(us,"cm") == 0) c = 1.0E6;
   if(strcmp(us,"A") == 0) c = 1.0E30;
   if(strcmp(us,"a.u.") == 0) c = 6.748321E30;
   if(c == 0.0) {
      fprintf(stderr,"unrecognised radius unit :%s\n",us);
      return(R);
   }
   if(dir < 0) c = 1/c;
   return(c * R);
}

double Tesla(E,us,dir)
double E;
char *us;
int dir;
{
double c;

   c = 0.0;
   if(strcmp(us,"T") == 0) c = 1.0;
   if(strcmp(us,"KG") == 0) c = 0.1;
   if(c == 0.0) {
      fprintf(stderr,"unrecognised Magnetic Field unit :%s\n",us);
      return(E);
   }
   if(dir < 0) c = 1/c;
   return(c * E);
}


double v(l,Lges,Sges)
double l,Lges,Sges;
{
double erg;

   erg = (4.0 * Sges - (2.0 * l + 1)) / ((2.0 * l - 1.0) * ( 2.0 * l + 3.0) * Sges * ( 2.0 * Lges - 1.0));
   return(erg);
}

double meJNJ(l,Lges,Sges,J)
double l,Lges,Sges,J;
{
double erg,JJ,LL,SS;

   JJ = J * (J + 1.0);
   LL = Lges * (Lges + 1.0);
   SS = Sges * (Sges + 1.0);
   if(J == 0) {
      erg = 1.0 + v(l,Lges,Sges) * 0.5 * Gs * (3.0 * JJ - 2.0 * LL - 3.0 * SS - 1.0);
   } else {
      erg = 3.0 * JJ - 2.0 * LL - 3.0 * SS - (3.0 * JJ - LL + SS) * (3.0 * JJ - LL - 3.0 * SS) / (4.0 * JJ);
      erg = erg * v(l,Lges,Sges) * 0.5 * Gs;
      erg = erg + (JJ + LL - SS)/(2.0 * JJ);
   }
   return(erg);
}

double sqr(x)
double x;
{
   return(x * x);
}

double meJAJ(Lges,Sges,J)
double Lges,Sges,J;
{
double erg,JJ,LL,SS;

   JJ = J * (J + 1.0);
   LL = Lges * (Lges + 1.0);
   SS = Sges * (Sges + 1.0);

   if(J == 0.0) return(Gs);
   erg = 1.0 + (Gs - 1) * (JJ + SS - LL) / (2.0 * JJ);
   return(erg);
}

double Bj(l,Lges,Sges,J)
double l,Lges,Sges,J;
{
double erg;

   erg = (sqr(Sges + Lges + 1.0) - J * J) * (J * J - sqr(Lges - Sges))/(6.0 * J);
   erg = erg * (1.0 + v(l,Lges,Sges) * (3 * J * J - Lges * (Lges + 1.0) - 3.0 * Sges * (Sges + 1.0)) / 2.0);
   return(erg);
}

double Calc_Bc(gj,J)
double gj,J;
{
return(-1.0 * B_core * (gj - 1) * J);
}

calc_ELSJ()
{
double Sg,Lg,x1,x2;
int i;

   Sg = 0.5 * ((2.0 * l + 1.0) - dabs(2.0 * l + 1.0 - Nel));
   Lg = Sg * dabs(2.0 * l + 1.0 - Nel);
   J_grund = Sg * dabs(2.0 * l - Nel);
   Fak_BN = 2.0 * MuB * Mu0 * R_mittel / (4.0 * Pi);
   for(i = 0; i < MAXTERM; i++) {
      Lges[i] = Lg;
      Sges[i] = Sg;
   }
   x1 = dabs(Lg - Sg);
   x2 = Lg + Sg;
   NTerm = (int) (x2 - x1 + 1.00001);
   for(i = 0; i < NTerm; i++) Jges[i] = (x1 + (double)i);
}


double Bi0(J)
double J;
{
double erg;
int n;

   n = 0;

   erg = Fak_BN * meJNJ(l,Lges[n],Sges[n],J) * J;
   erg = erg + Calc_Bc(meJAJ(Lges[n],Sges[n],J),J);
   return(erg);
}

double dexp(x)
double x;
{
double n,fak,px,erg,last;

   fak = 1.0; px = 1.0; n = 0.0; erg = 1.0; last = 0.0;
   while(erg != last) {
      last = erg;
      n = n + 1.0; fak = fak * n; px = px * x;
      erg = erg + px / fak;
      if(n > 100.0) break;
   }
   return(erg);
}

double Ehoch(x)
double x;
{
   if(x > 0.0) {
      if(x > 500.0) return(dexp(500.0)); else return(dexp(x));
   } else {
      if(x < -500.0) return(0.0); else return(dexp(x));
   }
}

/* initialize sum1, sum2, sum3 */
sum_init()
{
double J, gj, JBJ, JBJ1, f1, f2, Sg, Lg;
double dif;
int i;

   for(i = 0; i < MAXTERM; i++) {
      sum1[i] = 0.0;
      sum2[i] = 0.0;
      sum3[i] = 0.0;
   }

   for(i = 0; i < NTerm; i++) {

      J = Jges[i]; Lg = Lges[i]; Sg = Sges[i];

      gj = meJAJ(Lg,Sg,J);
      JBJ = Fak_BN * meJNJ(l,Lg,Sg,J);
      f1 = (2.0 * J + 1.0) * gj * MuB * (J + 1.0) * (JBJ * J + Calc_Bc(gj,J)) / (3.0 * k);
      sum1[i]= f1;

      if((i > 0) && (J != 0.0)) {
         JBJ1 = Fak_BN * Bj(l,Lg,Sg,J);
         dif = EJf[i] - EJf[i-1];
         f2 = MuB * JBJ1 / dif;
         sum2[i] = f2;
      }

      sum3[i] = (2.0 * J  + 1.0);
   }
}

double Beta(T)
double T;
{
double s1, s2, s3, kT, EXP_EJdkT;
int i;

   kT = k * T; s1 = 0.0; s2 = 0.0; s3 = 0.0;
   for(i = 0; i < NTerm; i++) {
      EXP_EJdkT = Ehoch(-1.0 * EJf[i] / kT);
/*
printf("i = %d , EJ/kt = %f , Exp(EJ/kT) = %f sum1[i] = %f\n",i,EJf[i]/kT,EXP_EJdkT,sum1[i]);
*/
      s1 = s1 + sum1[i] * EXP_EJdkT / T;
      if(i > 0) {
         s2 = s2 + sum2[i] * (Ehoch(-1.0 * EJf[i-1] / kT) - EXP_EJdkT);
/*
         printf("s2 = %-E , sum2[%d] = %-E \n",s2,i,sum2[i]);
*/
      }
      s3 = s3 + sum3[i] * EXP_EJdkT;
   }
/*
printf("s1 = %f , s2 = %f , s3 = %f\n",s1,s2,s3);
*/
   CurieTerm = 1.0E99;
   VanVleckTerm = 1.0E99;
   if(s3 != 0.0) {
      CurieTerm = s1 / s3;
      VanVleckTerm = s2 / s3;
      return(1.0 + (s1 - s2) / s3);
   } else {
      return(1.0E99);
   }
}

generate_spectrum(flag)
int flag;
{
int i;
double T;
float *curie, *vanvleck;

   if(flag == 2) {
      curie = (float *) calloc(_MAXSPCLEN,sizeof(float));
      vanvleck = (float *) calloc(_MAXSPCLEN,sizeof(float));
   }
   i = 0; T = Tbeg;
   while(T < Tend) {
      spc[i] = (float) Beta(T);
      err[i] = 0.0;
      tim[i] = (float) T; if(Flag_1dT == 1) tim[i] = (float) (1.0 / T);
      if(flag == 2) {
         curie[i] = CurieTerm;
         vanvleck[i] = VanVleckTerm;
      }
      T = T + Tstep; i = i + 1;
      if(i > _MAXSPCLEN) break;
   }
   unlink("Curie.spc");
   unlink("Curie.err");
   unlink("Curie.tim");
   unlink("VanVleck.spc");
   unlink("VanVleck.err");
   unlink("VanVleck.tim");
   writespec("BHF_out",spc,err,i,2,"B_HF-output");
   writespec("BHF_out.tim",tim,err,i,2,"B_HF-output");
   if(flag == 2) {
      writespec("Curie",curie,err,i,2,"Curie Term");
      writespec("Curie.tim",tim,err,i,2,"Curie Term");
      writespec("VanVleck",vanvleck,err,i,2,"Van Vleck Term");
      writespec("VanVleck.tim",tim,err,i,2,"Van Vleck Term");
      free(curie); free(vanvleck);
   }
}

plot_spectrum()
{
FILE *fp;

   SysCall("lise:ash bhf_out -v -pag");
   fp = fopen("Curie.spc","r");
   if(fp != NULL) {
      SysCall("lise:ash Curie -p 1 -1 -v -pag");
      fclose(fp);
   }
   fp = fopen("VanVleck.spc","r");
   if(fp != NULL) {
      SysCall("lise:ash VanVleck -p 1 -2 -v");
      fclose(fp);
   }
}


/* interactive Beta table */
table()
{
char *s;
int i;
double T,y;

   s = (char *) malloc(100);
   while(1 == 1) {
      printf("next temperature (or Just return to leave):"); fflush(stdout);
      fgets(s,80,stdin);
      if(s[0] == 10) break;
      sscanf(s,"%lf",&T);
      y = Beta(T);
      printf("CurieTerm = %-E , VanVleckTerm = %-E\n",CurieTerm,VanVleckTerm);
      printf("ß(%f) = %-E\n",T,y);
   }
   free(s);
}

show_faktor()
{
double J, gj, JBJ, JBJ1, f1, f2, Sg, Lg;
double dif;
int i;

   for(i = 0; i < MAXTERM; i++) {
      sum1[i] = 0.0;
      sum2[i] = 0.0;
      sum3[i] = 0.0;
   }

   for(i = 0; i < NTerm; i++) {

      J = Jges[i]; Lg = Lges[i]; Sg = Sges[i];

      gj = meJAJ(Lg,Sg,J);
      printf("gj(L=%2.2f,S=%2.2f,J=%2.2f) = %f\n",Lg,Sg,J,gj);
      JBJ = meJNJ(l,Lg,Sg,J);
      printf("<J|N|J>(l=%2.2f,L=%2.2f,S=%2.2f,J=%2.2f) = %f\n",l,Lg,Sg,J,JBJ);
      JBJ = Fak_BN * JBJ;
      f1 = (2.0 * J + 1.0) * gj * MuB * (J + 1.0) * (JBJ * J + Calc_Bc(gj,J)) / (3.0 * k);

      if((i > 0) && (J != 0.0)) {
         JBJ1 = Bj(l,Lg,Sg,J);
         printf("gj-1(l=%2.2f,L=%2.2f,S=%2.2f,J=%2.2f) = %f\n",l,Lg,Sg,J,JBJ1);
         JBJ1 = Fak_BN * JBJ1;
         dif = EJf[i] - EJf[i-1];
         f2 = MuB * JBJ1 / dif;
/*
printf("EJ[%d] = %-E , dif = %-E , JBJ1 = %-E\n",i,EJf[i],dif,JBJ1);
*/
      }
   }
}

double edit_double(text,x)
char *text;
double x;
{
char *s;
double y;

   s = (char *) malloc(100);
   printf("%s = %f ?",text,x); fflush(stdout);
   fgets(s,80,stdin);
   y = x;
   if(s[0] != 10) sscanf(s,"%lf",&y);
   free(s);
   return(y);
}

ask_data()
{
int i,n;
double x,y;
char s[4];

/*   gI = edit_double("Kern g Faktor",gI); */
   Nq = edit_double("Hauptquantenzahl",Nq);
   l = edit_double("Bahndrehimpuls (s=0,p=1,d=2,f=3)",l);
   Nel = edit_double("Anzahl der Elektronen in der Schale",Nel);
   x = ReCuR(R_mittel,"a.u.",-1);
   x = edit_double("<r-3> in atomaren Einheiten",x);
   R_mittel = ReCuR(x,"a.u.",1);
   x = Tesla(B_core,"T",-1);
   x = edit_double("Corepolarisation in Tesla",x);
   B_core = Tesla(x,"T",1);

   printf("calculating defaults ? (Y/N)"); fflush(stdout);
   fgets(s,2,stdin);
   if((s[0] == 'y') || (s[0] == 'Y')) calc_ELSJ();

   x = (double)NTerm;
   NTerm = (int) (0.00001 + edit_double("Anzahl der Terme",x));
   for(i = 0; i < NTerm; i++) {
      x = Energy(EJf[i],"K",-1);
      printf("%d : J = %f   S = %f   L = %f   EJ = %f\n",i,Jges[i],Sges[i],Lges[i],x);
      Jges[i] = edit_double("J",Jges[i]);
      Lges[i] = edit_double("L",Lges[i]);
      Sges[i] = edit_double("S",Sges[i]);
      x = edit_double("EJ in Kelvin",x);
      EJf[i] = Energy(x,"K",1);
   }
   printf("pre calculating sums...\n");
   sum_init();
   printf("ready\n");
}

calculate_data()
{
int i,n;
double x,y;
int flag;
char c,s[4];

   Tbeg = edit_double("Tbeg (Kelvin)",Tbeg);
   Tend = edit_double("Tend",Tend);
   Tstep = edit_double("Tstep",Tstep);
   flag = 1;
   printf("generate Curie curve and Van Vleck curve separately ? (Y/N)"); fflush(stdout);
   fgets(s,2,stdin); c = s[0];
   if((c == 'y') || (c == 'Y')) flag = 2;
   generate_spectrum(flag);
}

save_data()
{
FILE *fp;
int i,n;
double x;
char *s;

   s = (char *) malloc(100);
   printf("saving to >%s< ?",fnam); fflush(stdout);
   fgets(s,80,stdin);
   if(s[0] != 10) {s[strlen(s)-1] = 0; strcpy(fnam,s);}
   free(s);
   
   fp = fopen(fnam,"w");
   if(fp == NULL) {
      fprintf(stderr,"could not open file for write: %s\n",fnam);
      return();
   }

   fprintf(fp,"%f\n",gI);
   fprintf(fp,"%f\n",Nq);
   fprintf(fp,"%f\n",l);
   fprintf(fp,"%f\n",Nel);
   fprintf(fp,"%f\n",J_grund);
   fprintf(fp,"%f\n",R_mittel);
   fprintf(fp,"%f\n",B_core);
   fprintf(fp,"%d\n",NTerm);
   for(i = 0; i < MAXTERM; i++) {
      fprintf(fp,"%f\n",Jges[i]);
      fprintf(fp,"%f\n",Lges[i]);
      fprintf(fp,"%f\n",Sges[i]);
      x = Energy(EJf[i],"K",-1);
      fprintf(fp,"%f\n",x);
   }
   fclose(fp);
}

load_data()
{
int i,n;
char *s;

   s = (char *) malloc(100);
   printf("loading from >%s< ?",fnam); fflush(stdout);
   fgets(s,80,stdin);
   if(s[0] != 10) {s[strlen(s)-1] = 0; strcpy(fnam,s);}
   free(s);
   read_data();
}

read_data()
{
FILE *fp;
int i,n;
double x;

   fp = fopen(fnam,"r");
   if(fp == NULL) {
      fprintf(stderr,"could not open file for read: %s\n",fnam);
      printf("pre calculating sums...\n");
      sum_init();
      printf("ready\n");
      return();
   }
   fscanf(fp,"%lf",&gI);
   fscanf(fp,"%lf",&Nq);
   fscanf(fp,"%lf",&l);
   fscanf(fp,"%lf",&Nel);
   fscanf(fp,"%lf",&J_grund);
   fscanf(fp,"%lf",&R_mittel);
   fscanf(fp,"%lf",&B_core);
   fscanf(fp,"%d",&NTerm);
   for(i = 0; i < MAXTERM; i++) {
      fscanf(fp,"%lf",&Jges[i]);
      fscanf(fp,"%lf",&Lges[i]);
      fscanf(fp,"%lf",&Sges[i]);
      fscanf(fp,"%lf",&x);
      EJf[i] = Energy(x,"K",1);
   }
   fclose(fp);
   printf("pre calculating sums...\n");
   sum_init();
   printf("ready\n");
}

interactive()
{
int i, flag;
char *s,c;

   s = (char *) malloc(100);
   flag = 1;
   while(flag == 1) {
      printf("<S>ave (data)   <M>odify (data)   <C>alculate (spectrum)  <P>lot (spectrum)\n");
      printf("<L>oad (data)   <F>lag (1/T toggle)   <T>able (Beta(T))   <Q>uit\n");
      printf("<G> (gj,...)\n");
      fgets(s,80,stdin);
      c = s[0];
      switch(c) {
         case 'S': case 's': save_data(); break; 
         case 'L': case 'l': load_data(); break; 
         case 'M': case 'm': ask_data(); break;
         case 'C': case 'c': calculate_data(); break;
         case 'P': case 'p': plot_spectrum(); break;
         case 'F': case 'f': Flag_1dT = Flag_1dT ^ 1; break; 
         case 'T': case 't': table(); break;
         case 'Q': case 'q': flag = 0; break;
         case 'G': case 'g': show_faktor(); break;
         default: printf("unknown command >%s<\n",s);
      }
   }
   free(s);      
}

main(argc,argv)
int argc;
char *argv[];
{
int n,m,i,max;
char *z;

   spc = (float *) calloc(_MAXSPCLEN,sizeof(float));
   err = (float *) calloc(_MAXSPCLEN,sizeof(float));
   tim = (float *) calloc(_MAXSPCLEN,sizeof(float));
   EJf = (double *) calloc(MAXTERM,sizeof(double));
   Lges = (double *) calloc(MAXTERM,sizeof(double));
   Sges = (double *) calloc(MAXTERM,sizeof(double));
   Jges = (double *) calloc(MAXTERM,sizeof(double));
   sum1 = (double *) calloc(MAXTERM,sizeof(double));
   sum2 = (double *) calloc(MAXTERM,sizeof(double));
   sum3 = (double *) calloc(MAXTERM,sizeof(double));
   z = (char *) malloc(80);
   fnam = (char *) malloc(80);
   strcpy(fnam,"BHF_dat");

   for(i = 0; i < MAXTERM; i++) EJf[i] = (double) (i * 1000);

   R_mittel = ReCuR(6.796,"a.u.",1);
   EJf[1] = Energy(1600.0,"K",1);
   B_core = 10.0;
   calc_ELSJ();

   if(checkopt(argc,argv,"-r",z)) strcpy(fnam,z);
   if(checkopt(argc,argv,"-T",z)) Flag_1dT = 0;
   if(checkopt(argc,argv,"-Ts",z)) sscanf(z,"%lf",&Tbeg);
   if(checkopt(argc,argv,"-Te",z)) sscanf(z,"%lf",&Tend);
   if(checkopt(argc,argv,"-Ti",z)) sscanf(z,"%lf",&Tstep);

   read_data();
   if(checkopt(argc,argv,"-go",z)) {
      generate_spectrum(1);
   } else {
      help();
      interactive();
   }

   free(fnam); free(z);
   free(sum1); free(sum2); free(sum3);
   free(Jges); free(Sges); free(Lges);
   free(EJf); free(tim); free(err); free(spc);
}

SysCall(str)
char *str;
{
int i,n,m;

#ifdef AMIGA
struct FileHandle *input, *output;

input = (struct FileHandle *) Open("nil:",MODE_OLDFILE);
output = (struct FileHandle *) Open("nil:",MODE_NEWFILE);
   Execute(str,input,output);
Close(input);
Close(output);
#endif

#ifdef UNIX
   system(str);
#endif
}