برنامه تبدیل کارتزین به قطبی و قطبی به کارتزین در متلب

clc
clear all
disp ('for geodetic to cartesian press 1   for cartesino to zeodetic press 2')
ch=input  ('');
disp('select your elipsoid: WGS84=1   WGS72=2   GRS80=3')
el=input('');
 if el==1
    a=6378137;
    f=1/298.257223563;
end
if el==2
    a=6378135;
    f=1/298.26;
end
if el==3
    a=6378137;
    f=1/298.257222101;
end   
if el not 1:3 
    a=6378137;
    f=1/298.257223563;
end
if ch==1
disp('insert phi')
degp=input('degree=');
minp=input('minute=');
secp=input('second=');
secp=secp/3600;
minp=minp/60;
phie=degp+secp+minp;
phii=phie*(pi/180);
disp('insert landa')
degl=input('degree=');
minl=input('minute=');
secl=input('second=');
secl=secl/3600;
minl=minl/60;
phil=degl+secl+minl;
landa=phil*(pi/180);
disp('insert height')
hei=input('height=');
e=(2*f)-(f^2);
N=a/((1-(e*sin(phii)^2))^0.5);
X=(N+hei)*cos(phii)*cos(landa)
Y=(N+hei)*cos(phii)*sin(landa)
Z=(N*(1-e)+hei)*sin(phii)
else
X=input('Enter X=');
Y=input('Enter Y=');
Z=input('Enter Z=');
Landa=(atan(Y/X))*(180/pi);
disp('Landa=');
disp(Landa)
a=6378137;
f=1/298.257223563;
p=(((X^2)+(Y^2))^0.5);
e=((2*f)-(f^2));
phii(1)=atan((Z/p)*(1+(e/(1-e))));
N(1)=a/(1-(e*((sin(phii(1))^2)))^(1/2));
h(1)=(p/cos(phii(1)))-N(1);
phii(2)=atan((Z/p)*(1+(e*sin(phii(1))*(a/(sqrt(1-(e*(sin(phii(1)))^2))))/Z)));
N(2)=a/(1-(e*((sin(phii(2))^2)))^(1/2));
h(2)=(p/cos(phii(2)))-N(2);
for n=2:1000
 N(n)=a/(1-((e^2)*sin(phii(n-1)^2)))^(1/2);
 hi(n)=(p/cos(phii(n-1)))-N(n);
 phii(n)=atan((Z/p)*(1+(e*sin(phii(n-1)))*(a/(sqrt(1-(e*(sin(phii(n-1)))^2))))/Z));
   while (abs(phii(n)-phii(n-1))<=1/10^(1000000) & abs(N(n)-N(n-1))<=1/10^(1000000) & abs(hi(n)-hi(n-1))<=1/1000000 )
     break
   end
end
phii=phii(n)*(180/pi);
disp('phi=');disp(phii)
disp('height=');disp(hi(n))
end