function q0math
% Funtion for checking for q0 calculation math
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
	
	sep = repmat('_',1,80);
	
	syms h1 h2 positive % no loss of generality
	syms sin2t0 cos2t0 positive
	syms t real
 
	i1 = 1/(h1*cos(t)^2+h2*sin(t)^2);
	I1 = 4*limit(int(i1),'t',pi/2,'left');
	disp(sep), disp(['I1 = int(' char(i1) ',0,2*pi) = ']), pretty(I1)
	
	i2 = (cos2t0*cos(t)^2+sin2t0*sin(t)^2) / (h1*cos(t)^2+h2*sin(t)^2)^2;
	I2 = 4*limit(int(i2),'t',pi/2,'left');
	disp(sep), disp(['I2 = int(' char(i2) ',0,2*pi) = ']),	pretty(I2)
	
	i3 = ((h1-h2)^2*cos(t)^2*sin(t)^2) / (h1*cos(t)^2+h2*sin(t)^2)^2;
	I3 = 4*limit(int(i3),'t',pi/2,'left');
	disp(sep), disp(['I3 = int(' char(i3) ',0,2*pi) = ']),	pretty(I3)
	
	syms d2rF d2zF positive % no loss of generality
	syms drzF real
	J = [d2rF drzF ; drzF d2zF];
	[V,L] = eig(J);
	s = {h1 h2 sin2t0 cos2t0};
	r = {L(1,1) L(2,2) V(1,2)^2/(V(1,2)^2+V(2,2)^2) V(2,2)^2/(V(1,2)^2+V(2,2)^2)};
 
 N = h1*sin2t0+h2*cos2t0;
 subs(N,s,r)
	
	II1 = sqrt(simplify(subs(I1^2,s,r),'IgnoreAnalyticConstraints',true));
	disp(sep), disp('I1 ='),	pretty(II1)
	
	II2 = sqrt(simplify(subs(I2^2,s,r),'IgnoreAnalyticConstraints',true));
	disp(sep), disp('I2 ='),	pretty(II2)
	
	II3 = simplify(subs(I3,s,r),'IgnoreAnalyticConstraints',false);
	disp(sep), disp('I3 ='),	pretty(II3)
	
	I21 = subs(I2/I1,s,r);
	I21approx = simplify(taylor(I21,drzF,0,'order',3),'IgnoreAnalyticConstraints',true);
	disp(sep), disp('I2/I1 ='),	pretty(I2/I1)
	disp(' '), disp('to second order in drzF ='), pretty(I21approx)
	
	i1X = 1/(h1*cos(t)^2-h1*sin(t)^2);
	I1X = 4*limit(int(i1X),'t',pi/2,'left');
	disp(sep), disp(['I1X = int(' char(i1X) ',0,2*pi) = ']), pretty(I1X)
	
	i2X = (cos2t0*cos(t)^2+sin2t0*sin(t)^2) / (h1*cos(t)^2-h2*sin(t)^2)^2;
	I2X = 4*limit(int(i2X),'t',pi/2,'left');
	disp(sep), disp(['I2X = int(' char(i2X) ',0,2*pi) = ']),	pretty(I2X)
	
	i3 = ((h1-h2)^2*cos(t)^2*sin(t)^2) / (h1*cos(t)^2+h2*sin(t)^2)^2;
	I3 = 4*limit(int(i3),'t',pi/2,'left');
	disp(sep), disp(['I3 = int(' char(i3) ',0,2*pi) = ']),	pretty(I3)

 
 disp(sep), type q0math
end