/* FILE: uzeeman.cc            -*-Mode: c++-*-
 *
 * Uniform Zeeman (applied field) energy, derived from Oxs_Energy class.
 *
 */

#include <stdlib.h>
#include <vector>
#include "oc.h"
#include "nb.h"
#include "uzeeman.h"

// Oxs_Ext registration support
OXS_EXT_REGISTER(Oxs_UZeeman);

/* End includes */


// Constructor
Oxs_UZeeman::Oxs_UZeeman(
  const char* name,     // Child instance id
  Oxs_Director* newdtr, // App director
  const char* argstr)   // MIF input block parameters
  : Oxs_Energy(name,newdtr,argstr)
{
  // Process arguments
  // Get scale factor for all Hrange input
  REAL8m hmult = GetRealInitValue("multiplier",1.0);

  // Load Hrange data
  vector<string> hrange;
  FindRequiredInitValue("Hrange",hrange);

  // Build up Happ vector from Hrange data
  UINT4m step_count = 0;
  ThreeVector h0,h1;
  Oxs_SplitList range_list;
  for(UINT4m i=0;i<hrange.size();i++) {
    range_list.Split(hrange[i].c_str());
    if(range_list.Count() != 7) {
      char buf[1024];
      Oc_Snprintf(buf,sizeof(buf),
		  "Hrange[%d] value must be a 7 element list"
		  " (actual element count: %d)",
		  i,range_list.Count());
      throw Oxs_Ext::Error(this,buf);
    }
    BOOL err=0;
    h0.x = Nb_Atof(range_list[0],err);
    if(!err) h0.y = Nb_Atof(range_list[1],err);
    if(!err) h0.z = Nb_Atof(range_list[2],err);
    if(!err) h1.x = Nb_Atof(range_list[3],err);
    if(!err) h1.y = Nb_Atof(range_list[4],err);
    if(!err) h1.z = Nb_Atof(range_list[5],err);
    if(!err) {
      char* cptr;
      step_count = static_cast<UINT4m>(strtoul(range_list[6],&cptr,10));
      if(cptr == range_list[6] || *cptr!='\0') err=1;
    }
    if(err) {
      char buf[1024];
      Oc_Snprintf(buf,sizeof(buf),
		  "Hrange[%d] value is not a 7 element list"
		  " of 6 reals and an integer: %s",
		  i,hrange[i].c_str());
      throw Oxs_Ext::Error(this,buf);
    }
    // Note: There is a potential problem comparing Happ.back() to h0 if
    //  any floating point operations have occured on either.  If those
    //  operations are not exactly the same, there may be differences in
    //  round-off error.  Even if they are written the same in the source
    //  code, if the compiler is given license to re-order floating point
    //  ops, then the results may differ.  Also, there may be differences
    //  between values loaded from memory and those stored in registers.
    //  The code below attempts to work around these issues.
    if(step_count==0 || Happ.size()==0) {
      Happ.push_back(h0);
    } else {
      ThreeVector last_h1 = Happ.back();
      if(fabs(last_h1.x-h0.x)>fabs(h0.x)*REAL8_EPSILON*2
	 || fabs(last_h1.y-h0.y)>fabs(h0.y)*REAL8_EPSILON*2
	 || fabs(last_h1.z-h0.z)>fabs(h0.z)*REAL8_EPSILON*2) {
	Happ.push_back(h0);
      }
    }
    for(UINT4m j=1;j<step_count;j++) {
      REAL8m t = static_cast<REAL8m>(j)/static_cast<REAL8m>(step_count);
      Happ.push_back(((1-t)*h0) + (t*h1));
    }
    if(step_count>0) Happ.push_back(h1);
  }
  vector<ThreeVector>::iterator it = Happ.begin();
  while(it!=Happ.end()) {
    *it *= hmult;
    ++it;
  }
  DeleteInitValue("Hrange");
  VerifyAllInitArgsUsed();
}

Oxs_UZeeman::~Oxs_UZeeman()
{}

BOOL Oxs_UZeeman::Init()
{
  Bapp_output.Setup(this,InstanceName(),"B","mT",1,
		    &Oxs_UZeeman::Fill__Bapp_output,
		    director);
  Bappx_output.Setup(this,InstanceName(),"Bx","mT",1,
		     &Oxs_UZeeman::Fill__Bapp_output,
		     director);
  Bappy_output.Setup(this,InstanceName(),"By","mT",1,
		     &Oxs_UZeeman::Fill__Bapp_output,
		     director);
  Bappz_output.Setup(this,InstanceName(),"Bz","mT",1,
		     &Oxs_UZeeman::Fill__Bapp_output,
		     director);
  return Oxs_Energy::Init();
}

void Oxs_UZeeman::StageRequestCount
(unsigned int& min,
 unsigned int& max) const
{
  min = static_cast<unsigned int>(FieldCount());
  if(min<=1) max = UINT_MAX; // No upper limit
  else       max = min;      // Strict limit
}

ThreeVector Oxs_UZeeman::GetAppliedField(UINT4m stage_number) const
{
  ThreeVector H(0,0,0);
  UINT4m field_count = FieldCount();
  if(field_count>0) {
    if(stage_number>=field_count)  H=Happ[field_count-1];
    else                           H=Happ[stage_number];
  }

  return H;
}

void Oxs_UZeeman::GetEnergy
(const Oxs_SimState& state,
 Oxs_EnergyData& oed
 ) const
{
  UINT4m size = state.mesh->Size();
  if(size<1) return;

  // Use supplied buffer space, and reflect that use in oed.
  oed.energy = oed.energy_buffer;
  oed.field = oed.field_buffer;
  Oxs_MeshValue<REAL8m>& energy = *oed.energy_buffer;
  Oxs_MeshValue<ThreeVector>& field = *oed.field_buffer;

  ThreeVector H = GetAppliedField(state.stage_number);
  UINT4m i=0;
  do {
    field[i] = H;
    ++i;
  } while(i<size);

  ThreeVector vtemp(H);
  vtemp *= -MU0; // Constant so that
                /// vtemp * spin * Ms = cell energy density
  const Oxs_MeshValue<ThreeVector>& spin = state.spin;
  const Oxs_MeshValue<REAL8m>& Ms = *(state.Ms);
  i=0;
  do {
    energy[i] = Ms[i]*(vtemp*spin[i]);
    ++i;
  } while(i<size);

}

void
Oxs_UZeeman::Fill__Bapp_output(const Oxs_SimState& state)
{
  ThreeVector B = GetAppliedField(state.stage_number);

  B *= MU0*1000; // Report Bapp in mT

  if(Bapp_output.GetCacheRequestCount()>0) {
    Bapp_output.cache.state_id=0;
    Bapp_output.cache.value = sqrt(B.MagSq());
    Bapp_output.cache.state_id=state.Id();
  }

  if(Bappx_output.GetCacheRequestCount()>0) {
    Bappx_output.cache.state_id=0;
    Bappx_output.cache.value = B.x;
    Bappx_output.cache.state_id=state.Id();
  }

  if(Bappy_output.GetCacheRequestCount()>0) {
    Bappy_output.cache.state_id=0;
    Bappy_output.cache.value = B.y;
    Bappy_output.cache.state_id=state.Id();
  }

  if(Bappz_output.GetCacheRequestCount()>0) {
    Bappz_output.cache.state_id=0;
    Bappz_output.cache.value = B.z;
    Bappz_output.cache.state_id=state.Id();
  }

}