API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Dec 31st, 2014
527
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 20.07 KB | None | 0 0
raw download clone embed print report
  1. function [] = Earth_Orbit (TLE, t_step)
  2. %% INIT
  3. global mu RE
  4. mu = 398600.4418; %Mu km^3s^-2
  5. RE = 6378.1; %km
  6.  
  7. %% Inital State Determination
  8. %read in TLE
  9. [year,day,inc,RAAN,ecc,omega,M,n]=TLE_read(TLE);
  10.  
  11. %convert date to Julian
  12. JD=(juliandate(year,1,1,0,0,0))+day;
  13.  
  14. %Determine Inital Sate Vector
  15. COE_i=[JD,inc,RAAN,ecc,omega,M,n];
  16. [COE_i(8),COE_i(9),COE_i(10),COE_i(11)]=theta2Tah(COE_i(7), COE_i(4));
  17. [RVEC_i,VVEC_i]=COES2RV(COE_i(4),COE_i(2),COE_i(3),COE_i(5),COE_i(11),COE_i(9));
  18. [r_ECEF_i v_ECEF_i] = ECItoECEF(RVEC_i',VVEC_i',((JD-2451545)/36525),JD,-0.4399619,-0.140682,-0.33309,0,0,0);
  19. State = [RVEC_i(1),RVEC_i(2),RVEC_i(3),VVEC_i(1),VVEC_i(2),VVEC_i(3)];
  20.  
  21. %% Propagate Orbit
  22. options = odeset('RelTol',1e-8,'AbsTol',1e-8);
  23. [T,Y] = ode45(@CP,[0,t_step],State,options);
  24.  
  25. RVEC=Y(:,1:3);
  26. VVEC=Y(:,4:6);
  27.  
  28. rmag = norm(RVEC(length(RVEC),1:3));
  29. vmag = norm(VVEC(length(VVEC),1:3));
  30.  
  31. %% 3D plot
  32. earth_sphere
  33. hold on
  34. plot3(RVEC(:,1),RVEC(:,2),RVEC(:,3),'r')
  35. end
  36.  
  37. %% Orbit Propagator
  38. function [dr] = CP(t,r)
  39. rmag=norm(r(1:3));
  40.  
  41. %Define mu km^3s^-2
  42. global mu
  43.  
  44. dr(1,1)=r(4);
  45. dr(2,1)=r(5);
  46. dr(3,1)=r(6);
  47. dr(4,1)=(-(mu./rmag^3).*r(1));
  48. dr(5,1)=(-(mu./rmag^3).*r(2));
  49. dr(6,1)=(-(mu./rmag^3).*r(3));
  50. end
  51.  
  52. %%
  53.  
  54. function[year,day,inc,RAAN,ecc,omega,M,n]=TLE_read(TLE_file)
  55. TLE=fopen(TLE_file, 'r');
  56. line1=fgets(TLE);
  57. line2=fgets(TLE);
  58.  
  59. year=strcat('2','0',line1(19:20));
  60. year=str2double(year);
  61.  
  62. day=strcat(line1(21:32));
  63. day=str2double(day);
  64.  
  65. inc=strcat(line2(9:16));
  66. inc=str2double(inc);
  67.  
  68. RAAN=strcat(line2(18:25));
  69. RAAN=str2double(RAAN);
  70.  
  71. ecc=strcat('.',line2(27:33));
  72. ecc=str2double(ecc);
  73.  
  74. omega=strcat(line2(35:42));
  75. omega=str2double(omega);
  76.  
  77. M=strcat(line2(44:51));
  78. M=str2double(M);
  79.  
  80. n=strcat(line2(53:63));
  81. n=str2double(n);
  82.  
  83. end
  84.  
  85. function [T, a, h, theta] = theta2Tah(n, ecc)
  86.  
  87. global mu
  88.  
  89. T = 1/n;
  90. T = T*24*3600;
  91. a = (T*sqrt(mu)/2/pi)^(2/3);
  92. h = sqrt(a*(1 - ecc^2)*mu);
  93.  
  94. E=fzero(@(E) E-ecc*sin(E)-n,3);
  95. theta=2*atand(sqrt((1+ecc)/(1-ecc))*tan(E/2));
  96. end
  97.  
  98. function [RVEC, VVEC] = COES2RV(ecc, inc, RAAN, omega, theta, a)
  99.  
  100. global mu
  101.  
  102. %% RVEC VVEC PQW
  103. p=a*(1-ecc^2);
  104. RVEC_pqw = [(p*cosd(theta))/(1+ecc*cosd(theta));(p*sind(theta))/(1+ecc*cosd(theta));0];
  105. VVEC_pqw = [-sqrt(mu/p)*sind(theta);sqrt(mu/p)*(ecc+cosd(theta));0];
  106.  
  107. %% QxX Matrix (3x3)
  108. QXx = PeriToECIMatrix(RAAN, omega, inc);
  109.  
  110. %% RVEC_pqw VVEC_pqw => RVEC VVEC
  111. RVEC = (QXx * RVEC_pqw)';
  112. VVEC = (QXx * VVEC_pqw)';
  113.  
  114. end
  115.  
  116.  
  117. function QXx = PeriToECIMatrix(RAAN, omega, inc)
  118. QXx = [cosd(RAAN)*cosd(omega) - sind(RAAN)*sind(omega)*cosd(inc), ...
  119. - cosd(RAAN)*sind(omega) - sind(RAAN)*cosd(omega)*cosd(inc), ...
  120. sind(RAAN)*sind(inc); ...
  121. sind(RAAN)*cosd(omega) + cosd(RAAN)*sind(omega)*cosd(inc), ...
  122. -sind(RAAN)*sind(omega) + cosd(RAAN)*cosd(omega)*cosd(inc), ...
  123. -cosd(RAAN)*sind(inc);...
  124. sind(omega)*sind(inc), ...
  125. cosd(omega)*sind(inc), ...
  126. cosd(inc)];
  127. end
  128.  
  129. %% ECI to ECEF Function (Credit David Vallado)
  130. function [recef,vecef] = ECItoECEF ( reci,veci,ttt,jdut1,lod,xp,yp,eqeterms,ddpsi,ddeps);
  131. % function eci2ecef
  132. %
  133. % this function trsnforms a vector from the mean equator mean equniox frame
  134. % (j2000), to an earth fixed (ITRF) frame. the results take into account
  135. % the effects of precession, nutation, sidereal time, and polar motion.
  136. %
  137. % author : david vallado 719-573-2600 27 jun
  138. % 2002
  139.  
  140. %[prec,psia,wa,ea,xa] = precess ( ttt, '80' );
  141.  
  142. [deltapsi,trueeps,meaneps,omega,nut] = nutation(ttt,ddpsi,ddeps);
  143.  
  144. [st,stdot] = sidereal(jdut1,deltapsi,meaneps,omega,lod,eqeterms );
  145.  
  146. %[pm] = polarm(xp,yp,ttt,'80');
  147.  
  148.  
  149.  
  150. thetasa= 7.29211514670698e-05 * (1.0 - lod/86400.0 );
  151. omegaearth = [0; 0; thetasa;];
  152.  
  153. rpef = st'*reci;
  154. recef = rpef;
  155.  
  156.  
  157. vpef = st'*nut'*veci - cross( omegaearth,rpef );
  158. vecef = vpef;
  159. function [st,stdot] = sidereal (jdut1,deltapsi,meaneps,omega,lod,eqeterms );
  160. %
  161. % ----------------------------------------------------------------------------
  162. %
  163. % function sidereal
  164. %
  165. % this function calulates the transformation matrix that accounts for the
  166. % effects of sidereal time. Notice that deltaspi should not be moded to a
  167. % positive number because it is multiplied rather than used in a
  168. % trigonometric argument.
  169. %
  170. % author : david vallado 719-573-2600 25 jun
  171. % 2002-
  172. % ------------------------ find gmst --------------------------
  173. gmst= gstime( jdut1 );
  174.  
  175. % ------------------------ find mean ast ----------------------
  176. if (jdut1 > 2450449.5 ) & (eqeterms > 0)
  177. ast= gmst + deltapsi* cos(meaneps) ...
  178. + 0.00264*pi /(3600*180)*sin(omega) ...
  179. + 0.000063*pi /(3600*180)*sin(2.0 *omega);
  180. else
  181. ast= gmst + deltapsi* cos(meaneps);
  182. end
  183.  
  184. ast = rem (ast,2*pi);
  185. thetasa = 7.29211514670698e-05 * (1.0 - lod/86400.0 );
  186. omegaearth = thetasa;
  187.  
  188. %fprintf(1,'st gmst %11.7f ast %11.7f ome %11.7f \n',gmst*180/pi,ast*180/pi,omegaearth*180/pi );
  189.  
  190. st(1,1) = cos(ast);
  191. st(1,2) = -sin(ast);
  192. st(1,3) = 0.0;
  193. st(2,1) = sin(ast);
  194. st(2,2) = cos(ast);
  195. st(2,3) = 0.0;
  196. st(3,1) = 0.0;
  197. st(3,2) = 0.0;
  198. st(3,3) = 1.0;
  199.  
  200. % compute sidereal time rate matrix
  201. stdot(1,1) = -omegaearth * sin(ast);
  202. stdot(1,2) = -omegaearth * cos(ast);
  203. stdot(1,3) = 0.0;
  204. stdot(2,1) = omegaearth * cos(ast);
  205. stdot(2,2) = -omegaearth * sin(ast);
  206. stdot(2,3) = 0.0;
  207. stdot(3,1) = 0.0;
  208. stdot(3,2) = 0.0;
  209. stdot(3,3) = 0.0;
  210.  
  211. end
  212.  
  213. function [deltapsi, trueeps, meaneps, omega,nut] = nutation (ttt,ddpsi,ddeps );
  214.  
  215. deg2rad = pi/180.0;
  216.  
  217. [iar80,rar80] = iau80in; % coeff in deg
  218.  
  219. % ---- determine coefficients for iau 1980 nutation theory ----
  220. ttt2= ttt*ttt;
  221. ttt3= ttt2*ttt;
  222.  
  223. meaneps = -46.8150 *ttt - 0.00059 *ttt2 + 0.001813 *ttt3 + 84381.448;
  224. meaneps = rem( meaneps/3600.0 ,360.0 );
  225. meaneps = meaneps * deg2rad;
  226.  
  227. [ l, l1, f, d, omega, ...
  228. lonmer, lonven, lonear, lonmar, lonjup, lonsat, lonurn, lonnep, precrate ...
  229. ] = fundarg( ttt, '80' );
  230. %fprintf(1,'nut del arg %11.7f %11.7f %11.7f %11.7f %11.7f \n',l*180/pi,l1*180/pi,f*180/pi,d*180/pi,omega*180/pi );
  231.  
  232. deltapsi= 0.0;
  233. deltaeps= 0.0;
  234. for i= 106:-1: 1
  235. tempval= iar80(i,1)*l + iar80(i,2)*l1 + iar80(i,3)*f + ...
  236. iar80(i,4)*d + iar80(i,5)*omega;
  237. deltapsi= deltapsi + (rar80(i,1)+rar80(i,2)*ttt) * sin( tempval );
  238. deltaeps= deltaeps + (rar80(i,3)+rar80(i,4)*ttt) * cos( tempval );
  239. end
  240.  
  241. % --------------- find nutation parameters --------------------
  242. deltapsi = rem( deltapsi + ddpsi, 2.0 * pi );
  243. deltaeps = rem( deltaeps + ddeps, 2.0 * pi );
  244. trueeps = meaneps + deltaeps;
  245.  
  246. %fprintf(1,'meaneps %11.7f dp %11.7f de %11.7f te %11.7f ttt %11.7f \n',meaneps*180/pi,deltapsi*180/pi,deltaeps*180/pi,trueeps*180/pi, ttt );
  247.  
  248. cospsi = cos(deltapsi);
  249. sinpsi = sin(deltapsi);
  250. coseps = cos(meaneps);
  251. sineps = sin(meaneps);
  252. costrueeps = cos(trueeps);
  253. sintrueeps = sin(trueeps);
  254.  
  255. nut(1,1) = cospsi;
  256. nut(1,2) = costrueeps * sinpsi;
  257. nut(1,3) = sintrueeps * sinpsi;
  258. nut(2,1) = -coseps * sinpsi;
  259. nut(2,2) = costrueeps * coseps * cospsi + sintrueeps * sineps;
  260. nut(2,3) = sintrueeps * coseps * cospsi - sineps * costrueeps;
  261. nut(3,1) = -sineps * sinpsi;
  262. nut(3,2) = costrueeps * sineps * cospsi - sintrueeps * coseps;
  263. nut(3,3) = sintrueeps * sineps * cospsi + costrueeps * coseps;
  264.  
  265. % n1 = rot1mat( trueeps );
  266. % n2 = rot3mat( deltapsi );
  267. % n3 = rot1mat( -meaneps );
  268. % nut = n3*n2*n1
  269. end
  270.  
  271.  
  272.  
  273. function [iar80,rar80] = iau80in()
  274. % ----------------------------------------------------------------------------
  275. %
  276. % function iau80in
  277. %
  278. % this function initializes the nutation matricies needed for reduction
  279. % calculations. the routine needs the filename of the files as input.
  280. %
  281. % author : david vallado 719-573-2600 27 may
  282. % 2002
  283. % ------------------------ implementation -------------------
  284. % 0.0001" to rad
  285. convrt= 0.0001 * pi / (180*3600.0);
  286.  
  287. load nut80.dat;
  288.  
  289. iar80 = nut80(:,1:5);
  290. rar80 = nut80(:,6:9);
  291.  
  292. for i=1:106
  293. for j=1:4
  294. rar80(i,j)= rar80(i,j) * convrt;
  295. end
  296. end
  297. end
  298.  
  299.  
  300. function [ l, l1, f, d, omega, ...
  301. lonmer, lonven, lonear, lonmar, lonjup, lonsat, lonurn, lonnep, precrate ...
  302. ] = fundarg( ttt, opt )
  303. %
  304. % ----------------------------------------------------------------------------
  305. %
  306. % function fundarg
  307. %
  308. % this function calulates the delauany variables and planetary values for
  309. % several theories.
  310. %
  311. % author : david vallado 719-573-2600 16 jul 2004
  312.  
  313. iauhelp = 'n';
  314.  
  315. iaupnhelp = 'y';
  316.  
  317. show = 'n';
  318.  
  319. deg2rad = pi/180.0;
  320.  
  321. % ---- determine coefficients for iau 2000 nutation theory ----
  322. ttt2 = ttt*ttt;
  323. ttt3 = ttt2*ttt;
  324. ttt4 = ttt2*ttt2;
  325.  
  326. % ---- iau 2010 theory
  327. if opt == '10'
  328. % ------ form the delaunay fundamental arguments in deg
  329. l = 134.96340251 + ( 1717915923.2178 *ttt + ...
  330. 31.8792 *ttt2 + 0.051635 *ttt3 - 0.00024470 *ttt4 ) / 3600.0;
  331. l1 = 357.52910918 + ( 129596581.0481 *ttt - ...
  332. 0.5532 *ttt2 - 0.000136 *ttt3 - 0.00001149*ttt4 ) / 3600.0;
  333. f = 93.27209062 + ( 1739527262.8478 *ttt - ...
  334. 12.7512 *ttt2 + 0.001037 *ttt3 + 0.00000417*ttt4 ) / 3600.0;
  335. d = 297.85019547 + ( 1602961601.2090 *ttt - ...
  336. 6.3706 *ttt2 + 0.006593 *ttt3 - 0.00003169*ttt4 ) / 3600.0;
  337. omega= 125.04455501 + ( -6962890.5431 *ttt + ...
  338. 7.4722 *ttt2 + 0.007702 *ttt3 - 0.00005939*ttt4 ) / 3600.0;
  339.  
  340. % ------ form the planetary arguments in deg
  341. lonmer = 252.250905494 + 149472.6746358 *ttt;
  342. lonven = 181.979800853 + 58517.8156748 *ttt;
  343. lonear = 100.466448494 + 35999.3728521 *ttt;
  344. lonmar = 355.433274605 + 19140.299314 *ttt;
  345. lonjup = 34.351483900 + 3034.90567464 *ttt;
  346. lonsat = 50.0774713998 + 1222.11379404 *ttt;
  347. lonurn = 314.055005137 + 428.466998313*ttt;
  348. lonnep = 304.348665499 + 218.486200208*ttt;
  349. precrate= 1.39697137214*ttt + 0.0003086*ttt2;
  350. end;
  351.  
  352. % ---- iau 2000b theory
  353. if opt == '02'
  354. % ------ form the delaunay fundamental arguments in deg
  355. l = 134.96340251 + ( 1717915923.2178 *ttt ) / 3600.0;
  356. l1 = 357.52910918 + ( 129596581.0481 *ttt ) / 3600.0;
  357. f = 93.27209062 + ( 1739527262.8478 *ttt ) / 3600.0;
  358. d = 297.85019547 + ( 1602961601.2090 *ttt ) / 3600.0;
  359. omega= 125.04455501 + ( -6962890.5431 *ttt ) / 3600.0;
  360.  
  361. % ------ form the planetary arguments in deg
  362. lonmer = 0.0;
  363. lonven = 0.0;
  364. lonear = 0.0;
  365. lonmar = 0.0;
  366. lonjup = 0.0;
  367. lonsat = 0.0;
  368. lonurn = 0.0;
  369. lonnep = 0.0;
  370. precrate= 0.0;
  371. end;
  372.  
  373. % ---- iau 1996 theory
  374. if opt == '96'
  375. l = 134.96340251 + ( 1717915923.2178 *ttt + ...
  376. 31.8792 *ttt2 + 0.051635 *ttt3 - 0.00024470 *ttt4 ) / 3600.0;
  377. l1 = 357.52910918 + ( 129596581.0481 *ttt - ...
  378. 0.5532 *ttt2 - 0.000136 *ttt3 - 0.00001149*ttt4 ) / 3600.0;
  379. f = 93.27209062 + ( 1739527262.8478 *ttt - ...
  380. 12.7512 *ttt2 + 0.001037 *ttt3 + 0.00000417*ttt4 ) / 3600.0;
  381. d = 297.85019547 + ( 1602961601.2090 *ttt - ...
  382. 6.3706 *ttt2 + 0.006593 *ttt3 - 0.00003169*ttt4 ) / 3600.0;
  383. omega= 125.04455501 + ( -6962890.2665 *ttt + ...
  384. 7.4722 *ttt2 + 0.007702 *ttt3 - 0.00005939*ttt4 ) / 3600.0;
  385. % ------ form the planetary arguments in deg
  386. lonmer = 0.0;
  387. lonven = 181.979800853 + 58517.8156748 *ttt; % deg
  388. lonear = 100.466448494 + 35999.3728521 *ttt;
  389. lonmar = 355.433274605 + 19140.299314 *ttt;
  390. lonjup = 34.351483900 + 3034.90567464 *ttt;
  391. lonsat = 50.0774713998 + 1222.11379404 *ttt;
  392. lonurn = 0.0;
  393. lonnep = 0.0;
  394. precrate= 1.39697137214*ttt + 0.0003086*ttt2;
  395. end;
  396.  
  397. % ---- iau 1980 theory
  398. if opt == '80'
  399. l = 134.96298139 + ( 1717915922.6330 *ttt + ...
  400. 31.310 *ttt2 + 0.064 *ttt3 ) / 3600.0;
  401. l1 = 357.52772333 + ( 129596581.2240 *ttt - ...
  402. 0.577 *ttt2 - 0.012 *ttt3 ) / 3600.0;
  403. f = 93.27191028 + ( 1739527263.1370 *ttt - ...
  404. 13.257 *ttt2 + 0.011 *ttt3 ) / 3600.0;
  405. d = 297.85036306 + ( 1602961601.3280 *ttt - ...
  406. 6.891 *ttt2 + 0.019 *ttt3 ) / 3600.0;
  407. omega= 125.04452222 + ( -6962890.5390 *ttt + ...
  408. 7.455 *ttt2 + 0.008 *ttt3 ) / 3600.0;
  409. % ------ form the planetary arguments in deg
  410. lonmer = 252.3 + 149472.0 *ttt;
  411. lonven = 179.9 + 58517.8 *ttt; % deg
  412. lonear = 98.4 + 35999.4 *ttt;
  413. lonmar = 353.3 + 19140.3 *ttt;
  414. lonjup = 32.3 + 3034.9 *ttt;
  415. lonsat = 48.0 + 1222.1 *ttt;
  416. lonurn = 0.0;
  417. lonnep = 0.0;
  418. precrate= 0.0;
  419. end;
  420.  
  421. % ---- convert units to rad
  422. l = rem( l,360.0 ) * deg2rad; % rad
  423. l1 = rem( l1,360.0 ) * deg2rad;
  424. f = rem( f,360.0 ) * deg2rad;
  425. d = rem( d,360.0 ) * deg2rad;
  426. omega= rem( omega,360.0 ) * deg2rad;
  427.  
  428. lonmer= rem( lonmer,360.0 ) * deg2rad; % rad
  429. lonven= rem( lonven,360.0 ) * deg2rad;
  430. lonear= rem( lonear,360.0 ) * deg2rad;
  431. lonmar= rem( lonmar,360.0 ) * deg2rad;
  432. lonjup= rem( lonjup,360.0 ) * deg2rad;
  433. lonsat= rem( lonsat,360.0 ) * deg2rad;
  434. lonurn= rem( lonurn,360.0 ) * deg2rad;
  435. lonnep= rem( lonnep,360.0 ) * deg2rad;
  436. precrate= rem( precrate,360.0 ) * deg2rad;
  437. %iauhelp='y';
  438. if iauhelp == 'y'
  439. fprintf(1,'fa %11.7f %11.7f %11.7f %11.7f %11.7f deg \n',l*180/pi,l1*180/pi,f*180/pi,d*180/pi,omega*180/pi );
  440. fprintf(1,'fa %11.7f %11.7f %11.7f %11.7f deg \n',lonmer*180/pi,lonven*180/pi,lonear*180/pi,lonmar*180/pi );
  441. fprintf(1,'fa %11.7f %11.7f %11.7f %11.7f deg \n',lonjup*180/pi,lonsat*180/pi,lonurn*180/pi,lonnep*180/pi );
  442. fprintf(1,'fa %11.7f \n',precrate*180/pi );
  443. end;
  444. end
  445.  
  446. function gst = gstime(jdut1);
  447.  
  448. twopi = 2.0*pi;
  449. deg2rad = pi/180.0;
  450.  
  451. % ------------------------ implementation ------------------
  452. tut1= ( jdut1 - 2451545.0 ) / 36525.0;
  453.  
  454. temp = - 6.2e-6 * tut1 * tut1 * tut1 + 0.093104 * tut1 * tut1 ...
  455. + (876600.0 * 3600.0 + 8640184.812866) * tut1 + 67310.54841;
  456.  
  457. % 360/86400 = 1/240, to deg, to rad
  458. temp = rem( temp*deg2rad/240.0,twopi );
  459.  
  460. % ------------------------ check quadrants --------------------
  461. if ( temp < 0.0 )
  462. temp = temp + twopi;
  463. end
  464.  
  465. gst = temp;
  466. end
  467.  
  468. end
  469.  
  470. %% Earth Sphere (Credit Will Campbell)
  471. %http://www.mathworks.com/matlabcentral/fileexchange/27123-earth-sized-sphe
  472. %re-with-topography
  473.  
  474. function [xx,yy,zz] = earth_sphere(varargin)
  475. %EARTH_SPHERE Generate an earth-sized sphere.
  476. % [X,Y,Z] = EARTH_SPHERE(N) generates three (N+1)-by-(N+1)
  477. % matrices so that SURFACE(X,Y,Z) produces a sphere equal to
  478. % the radius of the earth in kilometers. The continents will be
  479. % displayed.
  480. %
  481. % [X,Y,Z] = EARTH_SPHERE uses N = 50.
  482. %
  483. % EARTH_SPHERE(N) and just EARTH_SPHERE graph the earth as a
  484. % SURFACE and do not return anything.
  485. %
  486. % EARTH_SPHERE(N,'mile') graphs the earth with miles as the unit rather
  487. % than kilometers. Other valid inputs are 'ft' 'm' 'nm' 'miles' and 'AU'
  488. % for feet, meters, nautical miles, miles, and astronomical units
  489. % respectively.
  490. %
  491. % EARTH_SPHERE(AX,...) plots into AX instead of GCA.
  492. %
  493. % Examples:
  494. % earth_sphere('nm') produces an earth-sized sphere in nautical miles
  495. %
  496. % earth_sphere(10,'AU') produces 10 point mesh of the Earth in
  497. % astronomical units
  498. %
  499. % h1 = gca;
  500. % earth_sphere(h1,'mile')
  501. % hold on
  502. % plot3(x,y,z)
  503. % produces the Earth in miles on axis h1 and plots a trajectory from
  504. % variables x, y, and z
  505.  
  506. % Clay M. Thompson 4-24-1991, CBM 8-21-92.
  507. % Will Campbell, 3-30-2010
  508. % Copyright 1984-2010 The MathWorks, Inc.
  509.  
  510. %% Input Handling
  511. [cax,args,nargs] = axescheck(varargin{:}); % Parse possible Axes input
  512. error(nargchk(0,2,nargs)); % Ensure there are a valid number of inputs
  513.  
  514. % Handle remaining inputs.
  515. % Should have 0 or 1 string input, 0 or 1 numeric input
  516. j = 0;
  517. k = 0;
  518. n = 50; % default value
  519. units = 'km'; % default value
  520. for i = 1:nargs
  521. if ischar(args{i})
  522. units = args{i};
  523. j = j+1;
  524. elseif isnumeric(args{i})
  525. n = args{i};
  526. k = k+1;
  527. end
  528. end
  529.  
  530. if j > 1 || k > 1
  531. error('Invalid input types')
  532. end
  533.  
  534. %% Calculations
  535.  
  536. % Scale factors
  537. Scale = {'km' 'm' 'mile' 'miles' 'nm' 'au' 'ft';
  538. 1 1000 0.621371192237334 0.621371192237334 0.539956803455724 6.6845871226706e-009 3280.839895};
  539.  
  540. % Identify which scale to use
  541. try
  542. myscale = 6378.1363*Scale{2,strcmpi(Scale(1,:),units)};
  543. catch %#ok<*CTCH>
  544. error('Invalid units requested. Please use m, km, ft, mile, miles, nm, or AU')
  545. end
  546.  
  547. % -pi <= theta <= pi is a row vector.
  548. % -pi/2 <= phi <= pi/2 is a column vector.
  549. theta = (-n:2:n)/n*pi;
  550. phi = (-n:2:n)'/n*pi/2;
  551. cosphi = cos(phi); cosphi(1) = 0; cosphi(n+1) = 0;
  552. sintheta = sin(theta); sintheta(1) = 0; sintheta(n+1) = 0;
  553.  
  554. x = myscale*cosphi*cos(theta);
  555. y = myscale*cosphi*sintheta;
  556. z = myscale*sin(phi)*ones(1,n+1);
  557.  
  558. %% Plotting
  559.  
  560. if nargout == 0
  561. cax = newplot(cax);
  562.  
  563. % Load and define topographic data
  564. load('topo.mat','topo','topomap1');
  565.  
  566. % Rotate data to be consistent with the Earth-Centered-Earth-Fixed
  567. % coordinate conventions. X axis goes through the prime meridian.
  568. % http://en.wikipedia.org/wiki/Geodetic_system#Earth_Centred_Earth_Fixed_.28ECEF_or_ECF.29_coordinates
  569. %
  570. % Note that if you plot orbit trajectories in the Earth-Centered-
  571. % Inertial, the orientation of the contintents will be misleading.
  572. topo2 = [topo(:,181:360) topo(:,1:180)]; %#ok<NODEF>
  573.  
  574. % Define surface settings
  575. props.FaceColor= 'texture';
  576. props.EdgeColor = 'none';
  577. props.FaceLighting = 'phong';
  578. props.Cdata = topo2;
  579.  
  580. % Create the sphere with Earth topography and adjust colormap
  581. surface(x,y,z,props,'parent',cax)
  582. colormap(topomap1)
  583.  
  584. % Replace the calls to surface and colormap with these lines if you do
  585. % not want the Earth's topography displayed.
  586. % surf(x,y,z,'parent',cax)
  587. % shading flat
  588. % colormap gray
  589.  
  590. % Refine figure
  591. axis equal
  592. xlabel(['X [' units ']'])
  593. ylabel(['Y [' units ']'])
  594. zlabel(['Z [' units ']'])
  595. view(127.5,30)
  596. else
  597. xx = x; yy = y; zz = z;
  598. end
  599. end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!