% MatLab R2015b
% Noncolaf.m file, v.21.04.2016

clearvars;
syms x y z
syms l1x l1y l1z l2x l2y l2z l3x l3y l3z
syms hx hy hz
syms h
syms lx ly lz mx my mz nx ny nz
syms xm ym zm x0 y0 z0 omega

%\\\\\\\\\\\ ENTER PARAMETERS ///////////

% Vector parameters l1 , l2 , l3 are parameterized by Euler angles in Radians:
% x = alpha angle [ 0 ; 2*pi ] , y = beta angle [ 0 ; pi ], z = gamma angle [ 0 ; 2*pi ]
% Projections of l1 , l2 , l3 on x , y , z axes:
% l1x , l1y , l1z
% l2x , l2y , l2z
% l3x , l3y , l3z

% Vector of external magnetic field is parameterized by hphi and htheta angles in Degrees:
% hphi [ 0 ; 360 deg ] , htheta [ - 90 deg ; 90 deg ]
% Projections of magnetic field vector H on x , y , z axes:
% Hx = H*cos(hphi)*cos(htheta) , Hy = H*sin(hphi)*cos(htheta) , Hz = H*sin(htheta)

% Magnetic field is expressed in kOe
% Frequency is expressed in GHz

% Hlow - lower boundary of magnetic field value
% Hhigh - higher boundary of magnetic field value
% delta - discontinuity of calculation for magnetic field value

% gamma - gyromagnetic ratio

% I1 , I2 , I3 - coefficients which are related to chi1 , chi2 , chi3 (static susceptibilities for magnetic field oriented
% parallel to vector parameters l1 , l2 , l3 respectively) in accordance with Andreev - Marchenko theory

% Uadd - relativistic potential energy of anisotropy which is expressed by
% l1x , l1y , l1z, l2x , l2y , l2z, l3x , l3y , l3z

% I1 , I2 , I3 , Uadd  have to be expressed in accordance with that
% magnetic field is in kOe and frequency is in GHz

hphi = 0; htheta = 0;
Hlow = 0;
Hhigh = 200;
delta = 1;

% test case of Mn3Al2Ge3O12
% gamma = 17.6;
% I1 = 1.42e-5; I2 = 1.42e-5; I3 = 7.99e-6;
% Uadd(l1x,l1y,l1z,l2x,l2y,l2z,l3x,l3y,l3z,h,hx,hy,hz) = 2/3^(1/2)*(l1x*l2x-l1y*l2y)+l1z^2-l2z^2;

% test case of LiCu2O2
% gamma = 17.6;
% I1 = 1.85e-7; I2 = 1.85e-7; I3 = 6.18e-8;
% Uadd(l1x,l1y,l1z,l2x,l2y,l2z,l3x,l3y,l3z,h,hx,hy,hz) = -1/2*l3z^2 - 99/200*l3y^2;

%test case of CsNiCl3
% gamma = 18.8;
% I1 = 8.77e-6; I2 = 8.77e-6; I3 = 9.75e-7;
% Uadd(l1x,l1y,l1z,l2x,l2y,l2z,l3x,l3y,l3z,h,hx,hy,hz) = 1/2*l3z^2;

 gamma = 17.6;
 I1 = 3e-6; I2 = 3e-6; I3 = 3e-6;
 Uadd(l1x,l1y,l1z,l2x,l2y,l2z,l3x,l3y,l3z,h,hx,hy,hz) = l3x^2+l3y^2;

%/////////// ENTER PARAMETERS \\\\\\\\\\\

%\\\\\\\\\\\ OUTPUT PARAMETERS ///////////

% All obtained data are saved in 3 files:

% 1) txt-file "Oscillation Eigenfrequencies.txt": data are arranged in 3 columns
%    H (kOe)   F1 (GHz)   F2 (GHz)   F3(GHz)

% H (kOe) - value of external magnetic field for which data are calculated in kOe
% F1 (GHz) - oscillation eigenfrequency 1 in GHz
% F2 (GHz) - oscillation eigenfrequency 2 in GHz
% F3 (GHz) - oscillation eigenfrequency 3 in GHz

% 2) txt-file "Static Properties.txt": data are arranged in 10 columns
%    H (kOe)   A1 (rad)   B1 (rad)   G1 (rad)   E (a.u.)   l1_paral   l2_paral   l3_paral   Chi_paral (a.u.)   Chi_perp (a.u.)

% A1 (rad) , B1 (rad) , G1 (rad) - alpha , beta , gamma angles for vector parameter l1 in Radians
% A2 = A1 , B2 = B1 , G2 = G1 + pi/2 - alpha , beta , gamma angles for vector parameter l2
% l3 = [ l1 , l2 ]
% E (a.u.) - value of full potential energy including magnetic and relativistic energy
% lj_paral - projection of lj vector on the field direction
% j = 1 , 2 , 3 - number of oscillation mode
% Chi_paral (a.u.) - static susceptibility among magnetic field vector
% Chi_perp (a.u.) - static susceptibility perpendicular to magnetic field vector

% 3) txt-file "Eigenvectors and Oscillating Magnetization Projections.txt": data are arranged in 25 columns
%    H (kOe)
%    Re[dA1]   Im[dA1]   Re[dB1]   Im[dB1]   Re[dG1]   Im[dG1]
%    Re[dA2]   Im[dA2]   Re[dB2]   Im[dB2]   Re[dG2]   Im[dG2]
%    Re[dA3]   Im[dA3]   Re[dB3]   Im[dB3]   Re[dG3]   Im[dG3]
%    M1_paral   M1_perp   M2_paral   M2_perp   M3_paral   M3_perp

% ( Re[dAj] + i*Im[dAj] , Re[dBj] + i*Im[dBj] , Re[dGj] + i*Im[dGj] ) - eigenvector for Fj - oscillation eigenfrequency
% j = 1 , 2 , 3 - number of oscillation mode
% Eigenvectors are normalized to unity

% Mj_paral^2 - average square of projection of oscillating magnetization on direction among magnetic field
% Mj_perp^2 - average square of projection of oscillating magnetization on plane perpendicular to magnetic field direction
% j = 1 , 2 , 3 - number of oscillation mode
% Mj_paral^2 + Mj_perp^2 = Mj^2 where Mj^2 is average square of oscillating magnetization
% Mj is normalized to unity

%/////////// OUTPUT PARAMETERS \\\\\\\\\\\

I1 = I1*gamma^2;
I2 = I2*gamma^2;
I3 = I3*gamma^2;

i = fix((Hhigh-Hlow)/delta);
hphi = hphi/180*pi; htheta = htheta/180*pi;
hxinit = cos(hphi)*cos(htheta); hyinit = sin(hphi)*cos(htheta); hzinit = sin(htheta);
arr(1:i+1,1:37) = 0;

fid1 = fopen('Static Properties.txt', 'w+');
fprintf(fid1, '# %s    %s    %s    %s    %s    %s    %s    %s    %s    %s\r\n', 'H (kOe)', 'A1 (rad)', 'B1 (rad)', 'G1 (rad)', 'E (a.u.)', 'l1_paral', 'l2_paral', 'l3_paral', 'Chi_paral (a.u.)', 'Chi_perp (a.u.)');
fid2 = fopen('Oscillation Eigenfrequencies.txt', 'w+');
fprintf(fid2, '# %s    %s    %s    %s\r\n', 'H (kOe)', 'F1 (GHz)', 'F2 (GHz)', 'F3 (GHz)');
fid3 = fopen('Eigenvectors and Oscillating Magnetization Projections.txt', 'w+');
fprintf(fid3, '# %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s    %s\r\n', 'H (kOe)', 'Re[dA1]', 'Im[dA1]', 'Re[dB1]', 'Im[dB1]', 'Re[dG1]', 'Im[dG1]', 'Re[dA2]', 'Im[dA2]', 'Re[dB2]', 'Im[dB2]', 'Re[dG2]', 'Im[dG2]', 'Re[dA3]', 'Im[dA3]', 'Re[dB3]', 'Im[dB3]', 'Re[dG3]', 'Im[dG3]', 'M1_paral (a.u.)', 'M1_perp (a.u.)', 'M2_paral (a.u.)', 'M2_perp (a.u.)', 'M3_paral (a.u.)', 'M3_perp (a.u.)');

hloc = Hlow-delta;

fx(x,y,z,lx,ly,lz) = cos(x)*cos(z)*lx-cos(x)*sin(z)*ly-sin(x)*cos(y)*sin(z)*lx-sin(x)*cos(y)*cos(z)*ly+sin(x)*sin(y)*lz;
fy(x,y,z,lx,ly,lz) = sin(x)*cos(z)*lx-sin(x)*sin(z)*ly+cos(x)*cos(y)*sin(z)*lx+cos(x)*cos(y)*cos(z)*ly-cos(x)*sin(y)*lz;
fz(x,y,z,lx,ly,lz) = sin(y)*sin(z)*lx+sin(y)*cos(z)*ly+cos(y)*lz;

gx(x,y,z,lx,ly,lz) = cos(x)*cos(z)*lx+sin(x)*cos(z)*ly-sin(x)*cos(y)*sin(z)*lx+cos(x)*cos(y)*sin(z)*ly+sin(y)*sin(z)*lz;
gy(x,y,z,lx,ly,lz) = -cos(x)*sin(z)*lx-sin(x)*sin(z)*ly-sin(x)*cos(y)*cos(z)*lx+cos(x)*cos(y)*cos(z)*ly+sin(y)*cos(z)*lz;
gz(x,y,z,lx,ly,lz) = sin(x)*sin(y)*lx-cos(x)*sin(y)*ly+cos(y)*lz;

l1xm(x,y,z) = cos(x)*cos(z)-sin(x)*cos(y)*sin(z);
l1ym(x,y,z) = sin(x)*cos(z)+cos(x)*cos(y)*sin(z);
l1zm(x,y,z) = sin(y)*sin(z);
l2xm(x,y,z) = -cos(x)*sin(z)-sin(x)*cos(y)*cos(z);
l2ym(x,y,z) = -sin(x)*sin(z)+cos(x)*cos(y)*cos(z);
l2zm(x,y,z) = sin(y)*cos(z);
l3xm(x,y,z) = sin(x)*sin(y);
l3ym(x,y,z) = -cos(x)*sin(y);
l3zm(x,y,z) = cos(y);

L1xm = matlabFunction(l1xm);
L1ym = matlabFunction(l1ym);
L1zm = matlabFunction(l1zm);
L2xm = matlabFunction(l2xm);
L2ym = matlabFunction(l2ym);
L2zm = matlabFunction(l2zm);
L3xm = matlabFunction(l3xm);
L3ym = matlabFunction(l3ym);
L3zm = matlabFunction(l3zm);

l1xm0 = L1xm(pi/6,pi/6,pi/6);
l1ym0 = L1ym(pi/6,pi/6,pi/6);
l1zm0 = L1zm(pi/6,pi/6,pi/6);
l2xm0 = L2xm(pi/6,pi/6,pi/6);
l2ym0 = L2ym(pi/6,pi/6,pi/6);
l2zm0 = L2zm(pi/6,pi/6,pi/6);
l3xm0 = L3xm(pi/6,pi/6,pi/6);
l3ym0 = L3ym(pi/6,pi/6,pi/6);
l3zm0 = L3zm(pi/6,pi/6,pi/6);

l1alpxm(x,y,z) = diff(l1xm(x,y,z),x);
l1betxm(x,y,z) = diff(l1xm(x,y,z),y);
l1gamxm(x,y,z) = diff(l1xm(x,y,z),z);
l1alpym(x,y,z) = diff(l1ym(x,y,z),x);
l1betym(x,y,z) = diff(l1ym(x,y,z),y);
l1gamym(x,y,z) = diff(l1ym(x,y,z),z);
l1alpzm(x,y,z) = diff(l1zm(x,y,z),x);
l1betzm(x,y,z) = diff(l1zm(x,y,z),y);
l1gamzm(x,y,z) = diff(l1zm(x,y,z),z);
l2alpxm(x,y,z) = diff(l2xm(x,y,z),x);
l2betxm(x,y,z) = diff(l2xm(x,y,z),y);
l2gamxm(x,y,z) = diff(l2xm(x,y,z),z);
l2alpym(x,y,z) = diff(l2ym(x,y,z),x);
l2betym(x,y,z) = diff(l2ym(x,y,z),y);
l2gamym(x,y,z) = diff(l2ym(x,y,z),z);
l2alpzm(x,y,z) = diff(l2zm(x,y,z),x);
l2betzm(x,y,z) = diff(l2zm(x,y,z),y);
l2gamzm(x,y,z) = diff(l2zm(x,y,z),z);
l3alpxm(x,y,z) = diff(l3xm(x,y,z),x);
l3betxm(x,y,z) = diff(l3xm(x,y,z),y);
l3gamxm(x,y,z) = diff(l3xm(x,y,z),z);
l3alpym(x,y,z) = diff(l3ym(x,y,z),x);
l3betym(x,y,z) = diff(l3ym(x,y,z),y);
l3gamym(x,y,z) = diff(l3ym(x,y,z),z);
l3alpzm(x,y,z) = diff(l3zm(x,y,z),x);
l3betzm(x,y,z) = diff(l3zm(x,y,z),y);
l3gamzm(x,y,z) = diff(l3zm(x,y,z),z);

scalar(lx,ly,lz,mx,my,mz) = lx*mx+ly*my+lz*mz;
Scalar = matlabFunction(scalar);
mixed(lx,ly,lz,mx,my,mz,nx,ny,nz) = (ly*mz-lz*my)*nx+(-lx*mz+lz*mx)*ny+(lx*my-ly*mx)*nz;

scalar1xx(x,y,z) = scalar(diff(l1xm(x,y,z),x),diff(l1ym(x,y,z),x),diff(l1zm(x,y,z),x),diff(l1xm(x,y,z),x),diff(l1ym(x,y,z),x),diff(l1zm(x,y,z),x));
scalar2xx(x,y,z) = scalar(diff(l2xm(x,y,z),x),diff(l2ym(x,y,z),x),diff(l2zm(x,y,z),x),diff(l2xm(x,y,z),x),diff(l2ym(x,y,z),x),diff(l2zm(x,y,z),x));
scalar3xx(x,y,z) = scalar(diff(l3xm(x,y,z),x),diff(l3ym(x,y,z),x),diff(l3zm(x,y,z),x),diff(l3xm(x,y,z),x),diff(l3ym(x,y,z),x),diff(l3zm(x,y,z),x));
scalar1yy(x,y,z) = scalar(diff(l1xm(x,y,z),y),diff(l1ym(x,y,z),y),diff(l1zm(x,y,z),y),diff(l1xm(x,y,z),y),diff(l1ym(x,y,z),y),diff(l1zm(x,y,z),y));
scalar2yy(x,y,z) = scalar(diff(l2xm(x,y,z),y),diff(l2ym(x,y,z),y),diff(l2zm(x,y,z),y),diff(l2xm(x,y,z),y),diff(l2ym(x,y,z),y),diff(l2zm(x,y,z),y));
scalar3yy(x,y,z) = scalar(diff(l3xm(x,y,z),y),diff(l3ym(x,y,z),y),diff(l3zm(x,y,z),y),diff(l3xm(x,y,z),y),diff(l3ym(x,y,z),y),diff(l3zm(x,y,z),y));
scalar1zz(x,y,z) = scalar(diff(l1xm(x,y,z),z),diff(l1ym(x,y,z),z),diff(l1zm(x,y,z),z),diff(l1xm(x,y,z),z),diff(l1ym(x,y,z),z),diff(l1zm(x,y,z),z));
scalar2zz(x,y,z) = scalar(diff(l2xm(x,y,z),z),diff(l2ym(x,y,z),z),diff(l2zm(x,y,z),z),diff(l2xm(x,y,z),z),diff(l2ym(x,y,z),z),diff(l2zm(x,y,z),z));
scalar3zz(x,y,z) = scalar(diff(l3xm(x,y,z),z),diff(l3ym(x,y,z),z),diff(l3zm(x,y,z),z),diff(l3xm(x,y,z),z),diff(l3ym(x,y,z),z),diff(l3zm(x,y,z),z));
scalar1xy(x,y,z) = scalar(diff(l1xm(x,y,z),x),diff(l1ym(x,y,z),x),diff(l1zm(x,y,z),x),diff(l1xm(x,y,z),y),diff(l1ym(x,y,z),y),diff(l1zm(x,y,z),y));
scalar2xy(x,y,z) = scalar(diff(l2xm(x,y,z),x),diff(l2ym(x,y,z),x),diff(l2zm(x,y,z),x),diff(l2xm(x,y,z),y),diff(l2ym(x,y,z),y),diff(l2zm(x,y,z),y));
scalar3xy(x,y,z) = scalar(diff(l3xm(x,y,z),x),diff(l3ym(x,y,z),x),diff(l3zm(x,y,z),x),diff(l3xm(x,y,z),y),diff(l3ym(x,y,z),y),diff(l3zm(x,y,z),y));
scalar1xz(x,y,z) = scalar(diff(l1xm(x,y,z),x),diff(l1ym(x,y,z),x),diff(l1zm(x,y,z),x),diff(l1xm(x,y,z),z),diff(l1ym(x,y,z),z),diff(l1zm(x,y,z),z));
scalar2xz(x,y,z) = scalar(diff(l2xm(x,y,z),x),diff(l2ym(x,y,z),x),diff(l2zm(x,y,z),x),diff(l2xm(x,y,z),z),diff(l2ym(x,y,z),z),diff(l2zm(x,y,z),z));
scalar3xz(x,y,z) = scalar(diff(l3xm(x,y,z),x),diff(l3ym(x,y,z),x),diff(l3zm(x,y,z),x),diff(l3xm(x,y,z),z),diff(l3ym(x,y,z),z),diff(l3zm(x,y,z),z));
scalar1yz(x,y,z) = scalar(diff(l1xm(x,y,z),y),diff(l1ym(x,y,z),y),diff(l1zm(x,y,z),y),diff(l1xm(x,y,z),z),diff(l1ym(x,y,z),z),diff(l1zm(x,y,z),z));
scalar2yz(x,y,z) = scalar(diff(l2xm(x,y,z),y),diff(l2ym(x,y,z),y),diff(l2zm(x,y,z),y),diff(l2xm(x,y,z),z),diff(l2ym(x,y,z),z),diff(l2zm(x,y,z),z));
scalar3yz(x,y,z) = scalar(diff(l3xm(x,y,z),y),diff(l3ym(x,y,z),y),diff(l3zm(x,y,z),y),diff(l3xm(x,y,z),z),diff(l3ym(x,y,z),z),diff(l3zm(x,y,z),z));

Scalar1xx = matlabFunction(scalar1xx);
Scalar2xx = matlabFunction(scalar2xx);
Scalar3xx = matlabFunction(scalar3xx);
Scalar1yy = matlabFunction(scalar1yy);
Scalar2yy = matlabFunction(scalar2yy);
Scalar3yy = matlabFunction(scalar3yy);
Scalar1zz = matlabFunction(scalar1zz);
Scalar2zz = matlabFunction(scalar2zz);
Scalar3zz = matlabFunction(scalar3zz);
Scalar1xy = matlabFunction(scalar1xy);
Scalar2xy = matlabFunction(scalar2xy);
Scalar3xy = matlabFunction(scalar3xy);
Scalar1xz = matlabFunction(scalar1xz);
Scalar2xz = matlabFunction(scalar2xz);
Scalar3xz = matlabFunction(scalar3xz);
Scalar1yz = matlabFunction(scalar1yz);
Scalar2yz = matlabFunction(scalar2yz);
Scalar3yz = matlabFunction(scalar3yz);

sc1xx = Scalar1xx(pi/6,pi/6,pi/6);
sc2xx = Scalar2xx(pi/6,pi/6,pi/6);
sc3xx = Scalar3xx(pi/6,pi/6,pi/6);
sc1yy = Scalar1yy(pi/6,pi/6,pi/6);
sc2yy = Scalar2yy(pi/6,pi/6,pi/6);
sc3yy = Scalar3yy(pi/6,pi/6,pi/6);
sc1zz = Scalar1zz(pi/6,pi/6,pi/6);
sc2zz = Scalar2zz(pi/6,pi/6,pi/6);
sc3zz = Scalar3zz(pi/6,pi/6,pi/6);
sc1xy = Scalar1xy(pi/6,pi/6,pi/6);
sc2xy = Scalar2xy(pi/6,pi/6,pi/6);
sc3xy = Scalar3xy(pi/6,pi/6,pi/6);
sc1xz = Scalar1xz(pi/6,pi/6,pi/6);
sc2xz = Scalar2xz(pi/6,pi/6,pi/6);
sc3xz = Scalar3xz(pi/6,pi/6,pi/6);
sc1yz = Scalar1yz(pi/6,pi/6,pi/6);
sc2yz = Scalar2yz(pi/6,pi/6,pi/6);
sc3yz = Scalar3yz(pi/6,pi/6,pi/6);

A = 1/gamma^2*(I1*sc1xx+I2*sc2xx+I3*sc3xx);
B = 1/gamma^2*(I1*sc1yy+I2*sc2yy+I3*sc3yy);
C = 1/gamma^2*(I1*sc1zz+I2*sc2zz+I3*sc3zz);
D = 1/gamma^2*(I1*sc1xy+I2*sc2xy+I3*sc3xy);
E = 1/gamma^2*(I1*sc1xz+I2*sc2xz+I3*sc3xz);
F = 1/gamma^2*(I1*sc1yz+I2*sc2yz+I3*sc3yz);

mixedxy(x,y,z,h,hx,hy,hz) = 2*h/gamma*(I1*mixed(diff(l1xm(x,y,z),x),diff(l1ym(x,y,z),x),diff(l1zm(x,y,z),x),diff(l1xm(x,y,z),y),diff(l1ym(x,y,z),y),diff(l1zm(x,y,z),y),hx,hy,hz)+I2*mixed(diff(l2xm(x,y,z),x),diff(l2ym(x,y,z),x),diff(l2zm(x,y,z),x),diff(l2xm(x,y,z),y),diff(l2ym(x,y,z),y),diff(l2zm(x,y,z),y),hx,hy,hz)+I3*mixed(diff(l3xm(x,y,z),x),diff(l3ym(x,y,z),x),diff(l3zm(x,y,z),x),diff(l3xm(x,y,z),y),diff(l3ym(x,y,z),y),diff(l3zm(x,y,z),y),hx,hy,hz));
mixedxz(x,y,z,h,hx,hy,hz) = 2*h/gamma*(I1*mixed(diff(l1xm(x,y,z),x),diff(l1ym(x,y,z),x),diff(l1zm(x,y,z),x),diff(l1xm(x,y,z),z),diff(l1ym(x,y,z),z),diff(l1zm(x,y,z),z),hx,hy,hz)+I2*mixed(diff(l2xm(x,y,z),x),diff(l2ym(x,y,z),x),diff(l2zm(x,y,z),x),diff(l2xm(x,y,z),z),diff(l2ym(x,y,z),z),diff(l2zm(x,y,z),z),hx,hy,hz)+I3*mixed(diff(l3xm(x,y,z),x),diff(l3ym(x,y,z),x),diff(l3zm(x,y,z),x),diff(l3xm(x,y,z),z),diff(l3ym(x,y,z),z),diff(l3zm(x,y,z),z),hx,hy,hz));
mixedyz(x,y,z,h,hx,hy,hz) = 2*h/gamma*(I1*mixed(diff(l1xm(x,y,z),y),diff(l1ym(x,y,z),y),diff(l1zm(x,y,z),y),diff(l1xm(x,y,z),z),diff(l1ym(x,y,z),z),diff(l1zm(x,y,z),z),hx,hy,hz)+I2*mixed(diff(l2xm(x,y,z),y),diff(l2ym(x,y,z),y),diff(l2zm(x,y,z),y),diff(l2xm(x,y,z),z),diff(l2ym(x,y,z),z),diff(l2zm(x,y,z),z),hx,hy,hz)+I3*mixed(diff(l3xm(x,y,z),y),diff(l3ym(x,y,z),y),diff(l3zm(x,y,z),y),diff(l3xm(x,y,z),z),diff(l3ym(x,y,z),z),diff(l3zm(x,y,z),z),hx,hy,hz));

Mixedxy = matlabFunction(mixedxy);
Mixedxz = matlabFunction(mixedxz);   %   Enter pi1/6 hloc hx00
Mixedyz = matlabFunction(mixedyz);

U0(l1x,l1y,l1z,l2x,l2y,l2z,l3x,l3y,l3z,h,hx,hy,hz) = -h^2/2*(I1+I2+I3)+h^2/2*(I1*(l1x*hx+l1y*hy+l1z*hz)^2+I2*(l2x*hx+l2y*hy+l2z*hz)^2+I3*(l3x*hx+l3y*hy+l3z*hz)^2);
U0add(x,y,z,h,hx,hy,hz) = U0(l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z),l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z),l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z),h,hx,hy,hz);
Uaddadd(x,y,z,h,hx,hy,hz) = Uadd(l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z),l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z),l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z),h,hx,hy,hz);
U1(x,y,z,h,hx,hy,hz) = U0add(x,y,z,h,hx,hy,hz) + Uaddadd(x,y,z,h,hx,hy,hz);  % Enter hloc hxinit

hx0(x,y,z,xm,ym,zm) = fx(x,y,z,gx(xm,ym,zm,hxinit,hyinit,hzinit),gy(xm,ym,zm,hxinit,hyinit,hzinit),gz(xm,ym,zm,hxinit,hyinit,hzinit));
hy0(x,y,z,xm,ym,zm) = fy(x,y,z,gx(xm,ym,zm,hxinit,hyinit,hzinit),gy(xm,ym,zm,hxinit,hyinit,hzinit),gz(xm,ym,zm,hxinit,hyinit,hzinit));
hz0(x,y,z,xm,ym,zm) = fz(x,y,z,gx(xm,ym,zm,hxinit,hyinit,hzinit),gy(xm,ym,zm,hxinit,hyinit,hzinit),gz(xm,ym,zm,hxinit,hyinit,hzinit));

Hx0 = matlabFunction(hx0);
Hy0 = matlabFunction(hy0);      % Enter pi/6, xmin
Hz0 = matlabFunction(hz0);

Uaddaddg(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = Uadd(fx(xm,ym,zm,gx(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z)),gy(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z)),gz(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z))),fy(xm,ym,zm,gx(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z)),gy(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z)),gz(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z))),fz(xm,ym,zm,gx(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z)),gy(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z)),gz(x0,y0,z0,l1xm(x,y,z),l1ym(x,y,z),l1zm(x,y,z))),fx(xm,ym,zm,gx(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z)),gy(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z)),gz(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z))),fy(xm,ym,zm,gx(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z)),gy(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z)),gz(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z))),fz(xm,ym,zm,gx(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z)),gy(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z)),gz(x0,y0,z0,l2xm(x,y,z),l2ym(x,y,z),l2zm(x,y,z))),fx(xm,ym,zm,gx(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z)),gy(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z)),gz(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z))),fy(xm,ym,zm,gx(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z)),gy(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z)),gz(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z))),fz(xm,ym,zm,gx(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z)),gy(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z)),gz(x0,y0,z0,l3xm(x,y,z),l3ym(x,y,z),l3zm(x,y,z))),h,hx,hy,hz);
U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = U0add(x,y,z,h,hx,hy,hz) + Uaddaddg(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz);

U2xx(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = diff(U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz),x,2);
U2yy(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = diff(U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz),y,2);
U2zz(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = diff(U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz),z,2);
U2xy(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = diff(U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz),x,y);
U2xz(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = diff(U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz),x,z);
U2yz(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz) = diff(U2(x,y,z,x0,y0,z0,xm,ym,zm,h,hx,hy,hz),y,z);

U2xx0 = matlabFunction(U2xx);
U2yy0 = matlabFunction(U2yy);
U2zz0 = matlabFunction(U2zz);
U2xy0 = matlabFunction(U2xy);       % Enter pi1/6 xmin hloc hx00
U2xz0 = matlabFunction(U2xz);
U2yz0 = matlabFunction(U2yz);

chixx(x,y,z) = I1*(1-(l1xm(x,y,z))^2)+I2*(1-(l2xm(x,y,z))^2)+I3*(1-(l3xm(x,y,z))^2);
chiyy(x,y,z) = I1*(1-(l1ym(x,y,z))^2)+I2*(1-(l2ym(x,y,z))^2)+I3*(1-(l3ym(x,y,z))^2);
chizz(x,y,z) = I1*(1-(l1zm(x,y,z))^2)+I2*(1-(l2zm(x,y,z))^2)+I3*(1-(l3zm(x,y,z))^2);
chixy(x,y,z) = -I1*l1xm(x,y,z)*l1ym(x,y,z)-I2*l2xm(x,y,z)*l2ym(x,y,z)-I3*l3xm(x,y,z)*l3ym(x,y,z);
chixz(x,y,z) = -I1*l1xm(x,y,z)*l1zm(x,y,z)-I2*l2xm(x,y,z)*l2zm(x,y,z)-I3*l3xm(x,y,z)*l3zm(x,y,z);
chiyz(x,y,z) = -I1*l1ym(x,y,z)*l1zm(x,y,z)-I2*l2ym(x,y,z)*l2zm(x,y,z)-I3*l3ym(x,y,z)*l3zm(x,y,z);

chiparal(x,y,z,hx,hy,hz) = scalar(chixx(x,y,z)*hx+chixy(x,y,z)*hy+chixz(x,y,z)*hz,chixy(x,y,z)*hx+chiyy(x,y,z)*hy+chiyz(x,y,z)*hz,chixz(x,y,z)*hx+chiyz(x,y,z)*hy+chizz(x,y,z)*hz,hx,hy,hz);
chiperp(x,y,z,hx,hy,hz) = (scalar(chixx(x,y,z)*hx+chixy(x,y,z)*hy+chixz(x,y,z)*hz,chixy(x,y,z)*hx+chiyy(x,y,z)*hy+chiyz(x,y,z)*hz,chixz(x,y,z)*hx+chiyz(x,y,z)*hy+chizz(x,y,z)*hz,chixx(x,y,z)*hx+chixy(x,y,z)*hy+chixz(x,y,z)*hz,chixy(x,y,z)*hx+chiyy(x,y,z)*hy+chiyz(x,y,z)*hz,chixz(x,y,z)*hx+chiyz(x,y,z)*hy+chizz(x,y,z)*hz)-(scalar(chixx(x,y,z)*hx+chixy(x,y,z)*hy+chixz(x,y,z)*hz,chixy(x,y,z)*hx+chiyy(x,y,z)*hy+chiyz(x,y,z)*hz,chixz(x,y,z)*hx+chiyz(x,y,z)*hy+chizz(x,y,z)*hz,hx,hy,hz))^2)^(1/2);
Chiparal = matlabFunction(chiparal);
Chiperp = matlabFunction(chiperp);

mx(x,y,z,xm,ym,zm,h,hx,hy,hz,omega) = 1i*gamma*omega*(I1*((l1alpym(xm,ym,zm)*l1zm(xm,ym,zm)-l1alpzm(xm,ym,zm)*l1ym(xm,ym,zm))*x+(l1betym(xm,ym,zm)*l1zm(xm,ym,zm)-l1betzm(xm,ym,zm)*l1ym(xm,ym,zm))*y+(l1gamym(xm,ym,zm)*l1zm(xm,ym,zm)-l1gamzm(xm,ym,zm)*l1ym(xm,ym,zm))*z)+I2*((l2alpym(xm,ym,zm)*l2zm(xm,ym,zm)-l2alpzm(xm,ym,zm)*l2ym(xm,ym,zm))*x+(l2betym(xm,ym,zm)*l2zm(xm,ym,zm)-l2betzm(xm,ym,zm)*l2ym(xm,ym,zm))*y+(l2gamym(xm,ym,zm)*l2zm(xm,ym,zm)-l2gamzm(xm,ym,zm)*l2ym(xm,ym,zm))*z)+I3*((l3alpym(xm,ym,zm)*l3zm(xm,ym,zm)-l3alpzm(xm,ym,zm)*l3ym(xm,ym,zm))*x+(l3betym(xm,ym,zm)*l3zm(xm,ym,zm)-l3betzm(xm,ym,zm)*l3ym(xm,ym,zm))*y+(l3gamym(xm,ym,zm)*l3zm(xm,ym,zm)-l3gamzm(xm,ym,zm)*l3ym(xm,ym,zm))*z))-gamma^2*h*(I1*(((l1alpxm(xm,ym,zm)*hx+l1alpym(xm,ym,zm)*hy+l1alpzm(xm,ym,zm)*hz)*l1xm(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1alpxm(xm,ym,zm))*x+((l1betxm(xm,ym,zm)*hx+l1betym(xm,ym,zm)*hy+l1betzm(xm,ym,zm)*hz)*l1xm(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1betxm(xm,ym,zm))*y+((l1gamxm(xm,ym,zm)*hx+l1gamym(xm,ym,zm)*hy+l1gamzm(xm,ym,zm)*hz)*l1xm(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1gamxm(xm,ym,zm))*z)+I2*(((l2alpxm(xm,ym,zm)*hx+l2alpym(xm,ym,zm)*hy+l2alpzm(xm,ym,zm)*hz)*l2xm(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2alpxm(xm,ym,zm))*x+((l2betxm(xm,ym,zm)*hx+l2betym(xm,ym,zm)*hy+l2betzm(xm,ym,zm)*hz)*l2xm(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2betxm(xm,ym,zm))*y+((l2gamxm(xm,ym,zm)*hx+l2gamym(xm,ym,zm)*hy+l2gamzm(xm,ym,zm)*hz)*l2xm(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2gamxm(xm,ym,zm))*z)+I3*(((l3alpxm(xm,ym,zm)*hx+l3alpym(xm,ym,zm)*hy+l3alpzm(xm,ym,zm)*hz)*l3xm(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3alpxm(xm,ym,zm))*x+((l3betxm(xm,ym,zm)*hx+l3betym(xm,ym,zm)*hy+l3betzm(xm,ym,zm)*hz)*l3xm(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3betxm(xm,ym,zm))*y+((l3gamxm(xm,ym,zm)*hx+l3gamym(xm,ym,zm)*hy+l3gamzm(xm,ym,zm)*hz)*l3xm(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3gamxm(xm,ym,zm))*z));
my(x,y,z,xm,ym,zm,h,hx,hy,hz,omega) = 1i*gamma*omega*(I1*((l1alpzm(xm,ym,zm)*l1xm(xm,ym,zm)-l1alpxm(xm,ym,zm)*l1zm(xm,ym,zm))*x+(l1betzm(xm,ym,zm)*l1xm(xm,ym,zm)-l1betxm(xm,ym,zm)*l1zm(xm,ym,zm))*y+(l1gamzm(xm,ym,zm)*l1xm(xm,ym,zm)-l1gamxm(xm,ym,zm)*l1zm(xm,ym,zm))*z)+I2*((l2alpzm(xm,ym,zm)*l2xm(xm,ym,zm)-l2alpxm(xm,ym,zm)*l2zm(xm,ym,zm))*x+(l2betzm(xm,ym,zm)*l2xm(xm,ym,zm)-l2betxm(xm,ym,zm)*l2zm(xm,ym,zm))*y+(l2gamzm(xm,ym,zm)*l2xm(xm,ym,zm)-l2gamxm(xm,ym,zm)*l2zm(xm,ym,zm))*z)+I3*((l3alpzm(xm,ym,zm)*l3xm(xm,ym,zm)-l3alpxm(xm,ym,zm)*l3zm(xm,ym,zm))*x+(l3betzm(xm,ym,zm)*l3xm(xm,ym,zm)-l3betxm(xm,ym,zm)*l3zm(xm,ym,zm))*y+(l3gamzm(xm,ym,zm)*l3xm(xm,ym,zm)-l3gamxm(xm,ym,zm)*l3zm(xm,ym,zm))*z))-gamma^2*h*(I1*(((l1alpxm(xm,ym,zm)*hx+l1alpym(xm,ym,zm)*hy+l1alpzm(xm,ym,zm)*hz)*l1ym(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1alpym(xm,ym,zm))*x+((l1betxm(xm,ym,zm)*hx+l1betym(xm,ym,zm)*hy+l1betzm(xm,ym,zm)*hz)*l1ym(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1betym(xm,ym,zm))*y+((l1gamxm(xm,ym,zm)*hx+l1gamym(xm,ym,zm)*hy+l1gamzm(xm,ym,zm)*hz)*l1ym(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1gamym(xm,ym,zm))*z)+I2*(((l2alpxm(xm,ym,zm)*hx+l2alpym(xm,ym,zm)*hy+l2alpzm(xm,ym,zm)*hz)*l2ym(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2alpym(xm,ym,zm))*x+((l2betxm(xm,ym,zm)*hx+l2betym(xm,ym,zm)*hy+l2betzm(xm,ym,zm)*hz)*l2ym(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2betym(xm,ym,zm))*y+((l2gamxm(xm,ym,zm)*hx+l2gamym(xm,ym,zm)*hy+l2gamzm(xm,ym,zm)*hz)*l2ym(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2gamym(xm,ym,zm))*z)+I3*(((l3alpxm(xm,ym,zm)*hx+l3alpym(xm,ym,zm)*hy+l3alpzm(xm,ym,zm)*hz)*l3ym(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3alpym(xm,ym,zm))*x+((l3betxm(xm,ym,zm)*hx+l3betym(xm,ym,zm)*hy+l3betzm(xm,ym,zm)*hz)*l3ym(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3betym(xm,ym,zm))*y+((l3gamxm(xm,ym,zm)*hx+l3gamym(xm,ym,zm)*hy+l3gamzm(xm,ym,zm)*hz)*l3ym(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3gamym(xm,ym,zm))*z));
mz(x,y,z,xm,ym,zm,h,hx,hy,hz,omega) = 1i*gamma*omega*(I1*((l1alpxm(xm,ym,zm)*l1ym(xm,ym,zm)-l1alpym(xm,ym,zm)*l1xm(xm,ym,zm))*x+(l1betxm(xm,ym,zm)*l1ym(xm,ym,zm)-l1betym(xm,ym,zm)*l1xm(xm,ym,zm))*y+(l1gamxm(xm,ym,zm)*l1ym(xm,ym,zm)-l1gamym(xm,ym,zm)*l1xm(xm,ym,zm))*z)+I2*((l2alpxm(xm,ym,zm)*l2ym(xm,ym,zm)-l2alpym(xm,ym,zm)*l2xm(xm,ym,zm))*x+(l2betxm(xm,ym,zm)*l2ym(xm,ym,zm)-l2betym(xm,ym,zm)*l2xm(xm,ym,zm))*y+(l2gamxm(xm,ym,zm)*l2ym(xm,ym,zm)-l2gamym(xm,ym,zm)*l2xm(xm,ym,zm))*z)+I3*((l3alpxm(xm,ym,zm)*l3ym(xm,ym,zm)-l3alpym(xm,ym,zm)*l3xm(xm,ym,zm))*x+(l3betxm(xm,ym,zm)*l3ym(xm,ym,zm)-l3betym(xm,ym,zm)*l3xm(xm,ym,zm))*y+(l3gamxm(xm,ym,zm)*l3ym(xm,ym,zm)-l3gamym(xm,ym,zm)*l3xm(xm,ym,zm))*z))-gamma^2*h*(I1*(((l1alpxm(xm,ym,zm)*hx+l1alpym(xm,ym,zm)*hy+l1alpzm(xm,ym,zm)*hz)*l1zm(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1alpzm(xm,ym,zm))*x+((l1betxm(xm,ym,zm)*hx+l1betym(xm,ym,zm)*hy+l1betzm(xm,ym,zm)*hz)*l1zm(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1betzm(xm,ym,zm))*y+((l1gamxm(xm,ym,zm)*hx+l1gamym(xm,ym,zm)*hy+l1gamzm(xm,ym,zm)*hz)*l1zm(xm,ym,zm)+(l1xm(xm,ym,zm)*hx+l1ym(xm,ym,zm)*hy+l1zm(xm,ym,zm)*hz)*l1gamzm(xm,ym,zm))*z)+I2*(((l2alpxm(xm,ym,zm)*hx+l2alpym(xm,ym,zm)*hy+l2alpzm(xm,ym,zm)*hz)*l2zm(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2alpzm(xm,ym,zm))*x+((l2betxm(xm,ym,zm)*hx+l2betym(xm,ym,zm)*hy+l2betzm(xm,ym,zm)*hz)*l2zm(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2betzm(xm,ym,zm))*y+((l2gamxm(xm,ym,zm)*hx+l2gamym(xm,ym,zm)*hy+l2gamzm(xm,ym,zm)*hz)*l2zm(xm,ym,zm)+(l2xm(xm,ym,zm)*hx+l2ym(xm,ym,zm)*hy+l2zm(xm,ym,zm)*hz)*l2gamzm(xm,ym,zm))*z)+I3*(((l3alpxm(xm,ym,zm)*hx+l3alpym(xm,ym,zm)*hy+l3alpzm(xm,ym,zm)*hz)*l3zm(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3alpzm(xm,ym,zm))*x+((l3betxm(xm,ym,zm)*hx+l3betym(xm,ym,zm)*hy+l3betzm(xm,ym,zm)*hz)*l3zm(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3betzm(xm,ym,zm))*y+((l3gamxm(xm,ym,zm)*hx+l3gamym(xm,ym,zm)*hy+l3gamzm(xm,ym,zm)*hz)*l3zm(xm,ym,zm)+(l3xm(xm,ym,zm)*hx+l3ym(xm,ym,zm)*hy+l3zm(xm,ym,zm)*hz)*l3gamzm(xm,ym,zm))*z));
Mx = matlabFunction(mx);
My = matlabFunction(my);
Mz = matlabFunction(mz);

mlong(x,y,z,xm,ym,zm,hx,hy,hz) = ((x*hx+y*hy+z*hz)^2+(xm*hx+ym*hy+zm*hz)^2)^(1/2)/(x^2+y^2+z^2+xm^2+ym^2+zm^2)^(1/2);
mtran(x,y,z,xm,ym,zm,hx,hy,hz) = (1-(mlong(x,y,z,xm,ym,zm,hx,hy,hz))^2)^(1/2);
Mlong = matlabFunction(mlong);
Mtran = matlabFunction(mtran);

for j = 1:i+1

hloc = hloc + delta;
arr(j,1) = hloc;

U10(x,y,z) = U1(x,y,z,hloc,hxinit,hyinit,hzinit);
U100 = matlabFunction(U10);
U1000 = @(w) U100(w(1), w(2), w(3));
gs = GlobalSearch;
problem = createOptimProblem('fmincon','x0',[pi,pi/2,pi],'objective',U1000,'lb',[0,0,0],'ub',[2*pi,pi,2*pi]);
[xmin,fmin,~,~,~] = run(gs,problem);

arr(j,2) = xmin(1); arr(j,3) = xmin(2); arr(j,4) = xmin(3);
arr(j,5) = fmin;

hx00 = Hx0(pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3));
hy00 = Hy0(pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3));
hz00 = Hz0(pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3));

arr(j,6) = Scalar(l1xm0,l1ym0,l1zm0,hx00,hy00,hz00);
arr(j,7) = Scalar(l2xm0,l2ym0,l2zm0,hx00,hy00,hz00);
arr(j,8) = Scalar(l3xm0,l3ym0,l3zm0,hx00,hy00,hz00);

arr(j,9) = Chiparal(pi/6,pi/6,pi/6,hx00,hy00,hz00);
arr(j,10) = Chiperp(pi/6,pi/6,pi/6,hx00,hy00,hz00);

d = Mixedxy(pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00);
e = Mixedxz(pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00);
f = Mixedyz(pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00);

U00gxx00 = U2xx0(pi/6,pi/6,pi/6,pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3),hloc,hx00,hy00,hz00);
U00gyy00 = U2yy0(pi/6,pi/6,pi/6,pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3),hloc,hx00,hy00,hz00);
U00gzz00 = U2zz0(pi/6,pi/6,pi/6,pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3),hloc,hx00,hy00,hz00);
U00gxy00 = U2xy0(pi/6,pi/6,pi/6,pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3),hloc,hx00,hy00,hz00);
U00gxz00 = U2xz0(pi/6,pi/6,pi/6,pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3),hloc,hx00,hy00,hz00);
U00gyz00 = U2yz0(pi/6,pi/6,pi/6,pi/6,pi/6,pi/6,xmin(1),xmin(2),xmin(3),hloc,hx00,hy00,hz00);

A1 = -A*B*C-2*D*E*F+A*F^2+B*E^2+C*D^2;
A0 = A*C*U00gyy00+B*C*U00gxx00+A*B*U00gzz00+2*D*E*U00gyz00+2*E*F*U00gxy00+2*D*F*U00gxz00-2*A*F*U00gyz00-F^2*U00gxx00-2*B*E*U00gxz00-E^2*U00gyy00-2*C*D*U00gxy00-D^2*U00gzz00+A*f^2+B*e^2+C*d^2+2*E*d*f-2*D*e*f-2*F*d*e;
B0 = -A*U00gyy00*U00gzz00-B*U00gxx00*U00gzz00-C*U00gxx00*U00gyy00-2*E*U00gxy00*U00gyz00-2*D*U00gxz00*U00gyz00-2*F*U00gxy00*U00gxz00+A*U00gyz00^2+2*F*U00gxx00*U00gyz00+B*U00gxz00^2+2*E*U00gyy00*U00gxz00+C*U00gxy00^2+2*D*U00gzz00*U00gxy00-U00gxx00*f^2-U00gyy00*e^2-U00gzz00*d^2+2*U00gxy00*e*f+2*U00gyz00*d*e-2*U00gxz00*d*f;
C0 = U00gxx00*U00gyy00*U00gzz00+2*U00gxy00*U00gxz00*U00gyz00-U00gxx00*U00gyz00^2-U00gyy00*U00gxz00^2-U00gzz00*U00gxy00^2;

Q = ((A0/A1)^2-3*B0/A1)/9;
R = (2*(A0/A1)^3-9*A0/A1*B0/A1+27*C0/A1)/54;
angle = 1/3*acos(R/Q^(3/2));
w1 = (-2*Q^(1/2)*cos(angle)-A0/A1/3)^(1/2);
w2 = (-2*Q^(1/2)*cos(angle+2*pi/3)-A0/A1/3)^(1/2);
w3 = (-2*Q^(1/2)*cos(angle-2*pi/3)-A0/A1/3)^(1/2);

arr(j,11) = w1/pi/2; arr(j,12) = w2/pi/2; arr(j,13) = w3/pi/2;

if w1/pi/2 < 1/10000
    
arr(j,32) = 0;
arr(j,33) = 0;

else
    
Am = U00gxx00 - w1^2*A;
Bm = U00gyy00 - w1^2*B;
Cm = U00gzz00 - w1^2*C;
Dm = U00gxy00 - w1^2*D;
Em = U00gxz00 - w1^2*E;
Fm = U00gyz00 - w1^2*F;

dm = w1*d;
em = w1*e;
fm = w1*f;

Matrix = [Am Dm+dm*1i Em+em*1i; Dm-dm*1i Bm Fm+fm*1i; Em-em*1i Fm-fm*1i Cm];
[Vmatr,Dmatr] = eig(Matrix);
flag = 1;
value = abs(Dmatr(1,1));

for s = 2:3
    if abs(Dmatr(s,s))<value
        flag = s;
        value = abs(Dmatr(s,s));
    end
end

arr(j,14) = real(Vmatr(1,flag)); arr(j,15) = imag(Vmatr(1,flag));
arr(j,16) = real(Vmatr(2,flag)); arr(j,17) = imag(Vmatr(2,flag));
arr(j,18) = real(Vmatr(3,flag)); arr(j,19) = imag(Vmatr(3,flag));

mx0 = Mx(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w1);
my0 = My(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w1);
mz0 = Mz(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w1);

arr(j,32) = Mlong(real(mx0),real(my0),real(mz0),imag(mx0),imag(my0),imag(mz0),hx00,hy00,hz00);
arr(j,33) = Mtran(real(mx0),real(my0),real(mz0),imag(mx0),imag(my0),imag(mz0),hx00,hy00,hz00);

end

if w2/pi/2 < 1/10000
    
arr(j,34) = 0;
arr(j,35) = 0;

else

Am = U00gxx00 - w2^2*A;
Bm = U00gyy00 - w2^2*B;
Cm = U00gzz00 - w2^2*C;
Dm = U00gxy00 - w2^2*D;
Em = U00gxz00 - w2^2*E;
Fm = U00gyz00 - w2^2*F;

dm = w2*d;
em = w2*e;
fm = w2*f;

Matrix = [Am Dm+dm*1i Em+em*1i; Dm-dm*1i Bm Fm+fm*1i; Em-em*1i Fm-fm*1i Cm];
[Vmatr,Dmatr] = eig(Matrix);
flag = 1;
value = abs(Dmatr(1,1));

for s = 2:3
    if abs(Dmatr(s,s))<value
        flag = s;
        value = abs(Dmatr(s,s));
    end
end

arr(j,20) = real(Vmatr(1,flag)); arr(j,21) = imag(Vmatr(1,flag));
arr(j,22) = real(Vmatr(2,flag)); arr(j,23) = imag(Vmatr(2,flag));
arr(j,24) = real(Vmatr(3,flag)); arr(j,25) = imag(Vmatr(3,flag));

mx0 = Mx(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w2);
my0 = My(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w2);
mz0 = Mz(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w2);

arr(j,34) = Mlong(real(mx0),real(my0),real(mz0),imag(mx0),imag(my0),imag(mz0),hx00,hy00,hz00);
arr(j,35) = Mtran(real(mx0),real(my0),real(mz0),imag(mx0),imag(my0),imag(mz0),hx00,hy00,hz00);

end

if w3/pi/2 < 1/10000
    
arr(j,36) = 0;
arr(j,37) = 0;

else

Am = U00gxx00 - w3^2*A;
Bm = U00gyy00 - w3^2*B;
Cm = U00gzz00 - w3^2*C;
Dm = U00gxy00 - w3^2*D;
Em = U00gxz00 - w3^2*E;
Fm = U00gyz00 - w3^2*F;

dm = w3*d;
em = w3*e;
fm = w3*f;

Matrix = [Am Dm+dm*1i Em+em*1i; Dm-dm*1i Bm Fm+fm*1i; Em-em*1i Fm-fm*1i Cm];
[Vmatr,Dmatr] = eig(Matrix);
flag = 1;
value = abs(Dmatr(1,1));

for s = 2:3
    if abs(Dmatr(s,s))<value
        flag = s;
        value = abs(Dmatr(s,s));
    end
end

arr(j,26) = real(Vmatr(1,flag)); arr(j,27) = imag(Vmatr(1,flag));
arr(j,28) = real(Vmatr(2,flag)); arr(j,29) = imag(Vmatr(2,flag));
arr(j,30) = real(Vmatr(3,flag)); arr(j,31) = imag(Vmatr(3,flag));

mx0 = Mx(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w3);
my0 = My(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w3);
mz0 = Mz(Vmatr(1,flag),Vmatr(2,flag),Vmatr(3,flag),pi/6,pi/6,pi/6,hloc,hx00,hy00,hz00,w3);

arr(j,36) = Mlong(real(mx0),real(my0),real(mz0),imag(mx0),imag(my0),imag(mz0),hx00,hy00,hz00);
arr(j,37) = Mtran(real(mx0),real(my0),real(mz0),imag(mx0),imag(my0),imag(mz0),hx00,hy00,hz00);

end

end

for j = 1:i+1

fprintf(fid1, '%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\r\n', arr(j,1:10));
fprintf(fid2, '%.7e\t%.7e\t%.7e\t%.7e\r\n', arr(j,1), arr(j,11), arr(j,12), arr(j,13));
fprintf(fid3, '%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\t%.7e\r\n', arr(j,1), arr(j,14), arr(j,15), arr(j,16), arr(j,17), arr(j,18), arr(j,19), arr(j,20), arr(j,21), arr(j,22), arr(j,23), arr(j,24), arr(j,25), arr(j,26), arr(j,27), arr(j,28), arr(j,29), arr(j,30), arr(j,31), arr(j,32), arr(j,33), arr(j,34), arr(j,35), arr(j,36), arr(j,37));

end

fclose(fid1);
fclose(fid2);
fclose(fid3);