برنامه تبدیل مختصات کارتزین به قطبی و بلعکس در متلب

شاید برای شما هم بعضی اوقات پیش اومده (مخصوصا دانشجویان رشته نقشه برداری) که مختصات ها را به هم تبدیل کنید

به روش های مختلفی میتوان مختصات کارتزین را به قطبی و قطبی را به کارتزین تبدیل کرد، که در این پست با استفاده از نرم افزار و زبان متلب سورسی را اماده کردیم که بتوان این کار را برای شما انجام دهد

این سورس هم به صورت عملی استفاده دارد و هم جهت یادگیری این زبان برنامه نویسی .

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