% Script version: 1.0, February 12, 2019
% (c) Vasiliy Glazkov, P.Kapitza Institute for Physical Problems

% calculating AFMR frequencies for collinear antiferromagnet
% at arbitrary field direction
% and at given anisotropy energy
% using hydrodynamic approach (Andreev&Marchenko, Sov.Phys.Usp 23,21 (1980))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  SETTING MODEL PARAMETERS  %%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%transverse susceptibility
chi=1.; %can be set to unity for AFMR calculations

%gyromagnetic ratio, GHz/kOe
gamma=1.;%for spin-only magnet (g=2) gamma=2.80 GHZ/kOe

%demonstartion sample of orthorhombic anisotropy
a1=4.;
a2=1.;

%parameters vector, to be used in energy calculations
%other parameters values to be added if necessary
parameters=[chi,gamma,a1,a2];

%field scan details
Hstart=0.;
Hstop=3.;
Hstep=0.01;

%field direction (polar coordinates)
Theta_H=0.*pi/180.;
Phi_H=90.*pi/180.;

%random minimum search tries
Ntries=30;

%save results
name="afmr.dat";

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% END OF MODEL PARAMETERS SETTINGS   %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%computing antiferromagnetic vector and its derivatives
%to avoid gimbal lock problem polar angles for l-vector 
%are counted from the rotated axes
function val=euler_matrix()
  %technical rotation of axes to avoid gimbal lock
  alpha=0.2*pi;
  beta=0.1*pi;
  gamma=0.15*pi;
  
  val=zeros(3,3);
  
  val(1,1)=cos(alpha)*cos(gamma)-sin(alpha)*cos(beta)*sin(gamma);
  val(2,1)=sin(alpha)*cos(gamma)+cos(alpha)*cos(beta)*sin(gamma);
  val(3,1)=sin(beta)*sin(gamma);
  val(1,2)=-cos(alpha)*sin(gamma)-sin(alpha)*cos(beta)*cos(gamma);
  val(2,2)=-sin(alpha)*sin(gamma)+cos(alpha)*cos(beta)*cos(gamma);
  val(3,2)=sin(beta)*cos(gamma);
  val(1,3)=sin(alpha)*sin(beta);
  val(2,3)=-cos(alpha)*sin(beta);
  val(3,3)=cos(beta);
endfunction 

function val=_l(angles)
  theta=angles(1);
  phi=angles(2);
  
  val=euler_matrix()*[sin(theta)*cos(phi);sin(theta)*sin(phi);cos(theta)];
endfunction

function val=_dldphi_i(angles,i)
  theta=angles(1);
  phi=angles(2);
  
  val=[0.;0.;0.];
  if(i==1)%polar angle derivative
    val=[cos(theta)*cos(phi);cos(theta)*sin(phi);-sin(theta)];
  endif
  if(i==2)%asimuthal angle derivative
    val=[-sin(theta)*sin(phi);sin(theta)*cos(phi);0];
  endif
  val=euler_matrix()*val;
endfunction

function val=_d2ldphi_ij(angles,i,j)
  theta=angles(1);
  phi=angles(2);
  val=[0.,0.,0.];
  if(i==1&&j==1)%double derivative over polar angle
    val=[-sin(theta)*cos(phi);-sin(theta)*sin(phi);-cos(theta)];
  endif  
 if(i==2&&j==2)%double derivative over asimuthal angle
    val=[-sin(theta)*cos(phi);-sin(theta)*sin(phi);0.];
  endif
 if((i==1&&j==2)||(i==2&&j==1))%mixed derivative
    val=[-cos(theta)*sin(phi);cos(theta)*cos(phi);0];
  endif
  val=euler_matrix()*val;
endfunction

% potential energy and Hessian matrix_type for single antiferromagnet
% assumed default anisotropy energy form
% U= (a1/2) l_x^2+(a2/2) l_y^2
% a1>a2>0, Z is the easy axis
%to be modified if additional anisotropies are included
function value=_U(parameters,vecH,angles)
  %preparing constants and parameters
  l=_l(angles);
  lx=l(1);
  ly=l(2);
  lz=l(3);  
  Hx=vecH(1);
  Hy=vecH(2);
  Hz=vecH(3);
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%% U(l) CALCULATIONS START   %%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
  chi=parameters(1);
  gamma=parameters(2);
  a1=parameters(3);
  a2=parameters(4);
              
  value=0.5*a1*(lx)^2+0.5*a2*(ly)^2;
  value+=0.5*chi*(dot(l,vecH))^2;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%% U(l) CALCULATIONS END     %%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  endfunction

function val=_d2U_dphi_ij(parameters,vecH,angles,i,j)
  %preparing constants and parameters
  
  l=_l(angles);
  dl_i=_dldphi_i(angles,i);
  dl_j=_dldphi_i(angles,j);
  d2l_ij=_d2ldphi_ij(angles,i,j);
 
  lx=l(1);
  ly=l(2);
  lz=l(3);
  Hx=vecH(1);
  Hy=vecH(2);
  Hz=vecH(3);
  dlx_i=dl_i(1);
  dly_i=dl_i(2);
  dlz_i=dl_i(3);
  dlx_j=dl_j(1);
  dly_j=dl_j(2);
  dlz_j=dl_j(3);
  d2lx_ij=d2l_ij(1);
  d2ly_ij=d2l_ij(2);
  d2lz_ij=d2l_ij(3);
  
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%% HESSIAN (ij) CALCULATIONS START   %%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  chi=parameters(1);
  gamma=parameters(2);
  a1=parameters(3);
  a2=parameters(4);  
  
  val=a1*(dlx_i*dlx_j+lx*d2lx_ij);
  val+=a2*(dly_i*dly_j+ly*d2ly_ij);
  val+=chi*(dot(dl_i,vecH)*dot(dl_j,vecH)+dot(l,vecH)*dot(d2l_ij,vecH));

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%   HESSIAN (ij) CALCULATIONS END   %%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

endfunction

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%                CALCULATIONS  ARE  BELOW                  %%%%%
%%%%%       NO CHANGES ARE REQUIRED HERE IN DEFAULT APPS       %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%recalculating parameters
%field direction
dirH=[sin(Theta_H)*cos(Phi_H);sin(Theta_H)*sin(Phi_H);cos(Theta_H)];

%initialising variables
H=Hstart;
out=[];
Harr=f1arr=f2arr=lxarr=lyarr=lzarr=[];

angles_min=[0.1,0.1];

while(H<Hstop)
  %magnetic field vector
  vecH=H*dirH;
  
  %looking for equilibrium position
  %Search over random grid for best approximation
  Umin=_U(parameters,vecH,angles_min);

  for(n=1:1:Ntries)
    l=[normrnd(0,1);normrnd(0,1);normrnd(0,1)];
    l=l/sqrt(l(1)^2+l(2)^2+l(3)^2);
    angles(1)=acos(l(3));
    angles(2)=0.;
    if(abs(sin(angles(1)))>0.01)%not close to Z axis
      _cos=l(1)/sin(angles(1));
      _sin=l(2)/sin(angles(1));
      phi=acos(_cos);
      if(_sin<0)
        phi=-phi;
      endif
    endif  
    angles(2)=phi;
    
    U=_U(parameters,vecH,angles);
    
    if(U<Umin)
      angles_min=angles;
      Umin=U;
    endif
  endfor
  
  %search for minimum accurately
  angles_min=fminsearch(@(angles) _U(parameters,vecH,angles),angles_min);
  l=_l(angles_min);
  
  lxarr=[lxarr,l(1)];
  lyarr=[lyarr,l(2)];
  lzarr=[lzarr,l(3)];
  
  %solving for eigenfrequencies
  %Af^4-Bf^2+C=0
  
  dl1=_dldphi_i(angles_min,1);
  dl2=_dldphi_i(angles_min,2);
  dU1=_d2U_dphi_ij(parameters,vecH,angles_min,1,1);
  dU2=_d2U_dphi_ij(parameters,vecH,angles_min,2,2);
  dU12=_d2U_dphi_ij(parameters,vecH,angles_min,1,2);
  
  A=(chi/gamma^2)^2*dot(dl1,dl1)*dot(dl2,dl2);
  A+=(chi/gamma^2)^2*(dot(dl1,dl2))^2;  
  B=(chi/gamma^2)*dot(dl1,dl1)*dU2+(chi/gamma^2)*dot(dl2,dl2)*dU1;
  B+=4*(chi/gamma)^2*(dot(dl1,cross(dl2,vecH)))^2;
  B-=2*(chi/gamma^2)*dot(dl1,dl2)*dU12;
  C=dU1*dU2-(dU12)^2;
  
  D=B^2-4*A*C;
  
  f1=f2=NaN;
  if(D>=0)
    sqf1=(B+sqrt(D))/(2*A);
    sqf2=(B-sqrt(D))/(2*A);
    if(sqf1>=0)
      f1=sqrt(sqf1);
    endif
    if(sqf2>=0)
      f2=sqrt(sqf2);
    endif
  endif

  f1arr=[f1arr,f1];
  f2arr=[f2arr,f2];
  Harr=[Harr,H];
  
  plot(Harr,f1arr,Harr,f2arr);
  refresh;
  
  H+=Hstep;
endwhile 

out=transpose([Harr;f1arr;f2arr;lxarr;lyarr;lzarr]);
save("-text",name,"out");