% 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))