API tools faq
paste
Advertisement
Guest User

Matlab/Octave/Compose

a guest
Mar 17th, 2019
99
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Octave 12.69 KB | None | 0 0
raw download clone embed print report
  1. clear %Because when variables, such as a, b, ti[me] are changed, sizes of arrays get changed, and it's inconvenient to stumble over "array of wrong size" all over again
  2. %When translating from Maple to Matlab/Octave/OML, replace all := with =, and then replace all : with ;
  3. REa=6.3781;REb=6.3568;RM=1.73814;RS=696;%Constants
  4. ESa=152098.232;ESb=147098.290;ESan=0;EStilt=pi*(23+26/60+21/3600)/180;%Pi in Maple, case matters
  5. ST=8766.1536;MT=27.3*24;SMtilt=pi*5.1/180;
  6. EMa=405.696;EMb=363.104;
  7. %a='a';b='b'; in Maple, equations are symbolic, and variables such as a, b
  8. %and time need to be given values only at the end
  9. a=(0:0.05:2)*pi;b=0.0:0.02:1.0;ti=-2.5:0.05:2.5;tim=1:length(ti);%for "array index"%b should be between 0 and 1, but poles are difficult to evaluate, division by 0
  10. tir=repmat (ti', [1 length(a) length(b)]);
  11. %Sun movement relative to Earth
  12. SCX(tim)=ESa*cos(ESan+2*pi*ti/ST)*cos(EStilt);%Center of the Sun
  13. SX=zeros(length(ti),length(a),length(b));%Surface of the Sun
  14. for index=1:length(ti)%bad fix, but 3D arrays in Matlab get irritating
  15.     SX(index,:,:)=RS*cos(a')*cos(b*2*pi)+SCX(index);
  16. end
  17. SCY(tim)=ESb*sin(ESan+((2.*pi.*ti)/ST));
  18. %SY(tim)=RS*cos(a')*sin(b*2*pi)+SCY(tim);
  19. for index=1:length(ti)
  20.     SY(index,:,:)=RS*cos(a')*sin(b*2*pi)+SCY(index);
  21. end
  22. SCZ(tim)=ESa*cos(ESan+((2*pi*ti)/ST))*sin(EStilt); %???
  23. %SZ(tim)=RS*sin(a)+SCZ(tim);
  24. for index=1:length(ti)
  25.     SZ(index,:,:)=RS*sin(a')*ones(1,length(b))+SCZ(index);
  26. end
  27.  
  28. try
  29.     SXd(tim,:,:)=SX(tim,:,:).*cos(pi*ti'/12)+SY(tim,:,:).*sin(pi*ti'/12); %Matlab
  30.     SYd(tim,:,:)=SY(tim,:,:).*cos(pi*ti'/12)-SX(tim,:,:).*sin(pi*ti'/12);
  31. catch
  32.     SXd(tim,:,:)=SX(tim,:,:).*cos(pi*tir/12)+SY(tim,:,:).*sin(pi*tir/12); %Compose
  33.     SYd(tim,:,:)=SY(tim,:,:).*cos(pi*tir/12)-SX(tim,:,:).*sin(pi*tir/12);
  34. end
  35.  
  36. SCXd(tim)=SCX(tim).*cos(pi*ti/12)+SCY(tim).*sin(pi*ti/12);
  37. SCYd(tim)=SCY(tim).*cos(pi*ti/12)-SCX(tim).*sin(pi*ti/12);
  38.  
  39. %Moon movement relative to Earth
  40. MXel(tim)=EMb*cos(2*pi*ti/MT); MYel(tim)=EMa*sin(2*pi*ti/MT);
  41. MXSM(tim)=MXel(tim); MYSM(tim)=MYel(tim)*cos(SMtilt); MZSM(tim)=MYel(tim)*sin(SMtilt);
  42. MXSME(tim)=MXSM(tim)*cos(ESan)+MYSM(tim)*sin(ESan); MYSME(tim)=MYSM(tim)*cos(ESan)-MXSM(tim)*sin(ESan); MZSME(tim)=MZSM(tim);
  43. MCX(tim)=MXSME(tim)*cos(EStilt)-MZSME(tim)*sin(EStilt); %Center of the Moon
  44. %MX(tim)=RM.*cos(a')*cos(b*2*pi)+MCX(tim);
  45. for index=1:length(ti)
  46.     MX(index,:,:)=RM*cos(a')*cos(b*2*pi)+MCX(index);
  47. end
  48. MCY(tim)=MYSME(tim);
  49. %MY(tim)=RM.*cos(a).*sin(b.*2.*pi)+MCY(tim);
  50. for index=1:length(ti)
  51.     MY(index,:,:)=RM*cos(a')*sin(b*2*pi)+MCY(index);
  52. end
  53. MCZ(tim)=MZSME(tim).*cos(EStilt)+MXSME(tim).*sin(EStilt);
  54. %MZ(tim)=RM.*sin(a)+MCZ(tim);
  55. for index=1:length(ti)
  56.     MZ(index,:,:)=RM*sin(a')*ones(1,length(b))+MCZ(index);
  57. end
  58.  
  59. try
  60.     MXd(tim,:,:)=MX(tim,:,:).*cos(pi*ti'/12)+MY(tim,:,:).*sin(pi*ti'/12); %Matlab
  61.     MYd(tim,:,:)=MY(tim,:,:).*cos(pi*ti'/12)-MX(tim,:,:).*sin(pi*ti'/12);
  62. catch
  63.     MXd(tim,:,:)=MX(tim,:,:).*cos(pi*tir/12)+MY(tim,:,:).*sin(pi*tir/12); %Compose
  64.     MYd(tim,:,:)=MY(tim,:,:).*cos(pi*tir/12)-MX(tim,:,:).*sin(pi*tir/12);
  65. end
  66. MCXd(tim)=MCX(tim).*cos(pi*ti/12)+MCY(tim).*sin(pi*ti/12);
  67. MCYd(tim)=MCY(tim).*cos(pi*ti/12)-MCX(tim).*sin(pi*ti/12);
  68.  
  69. %Distance between the Moon and the Sun
  70. MSx(tim)=MCXd(tim)-SCXd(tim);
  71. MSy(tim)=MCYd(tim)-SCYd(tim);
  72. MSz(tim)=MCZ(tim)-SCZ(tim);
  73. MS(tim)=((MSx(tim)).^2+(MSy(tim)).^2+(MSz(tim)).^2).^(0.5);
  74. %MSa(tim)=arcsin((RS-RM)./MS(tim)); MSo(tim)=arcsin((RS+RM)./MS(tim));%atan
  75. %too instead of arctan
  76. MSa(tim)=asin((RS-RM)./MS(tim)); MSo(tim)=asin((RS+RM)./MS(tim));
  77.  
  78. %Moon shadow on Earth
  79. px1(tim)=0;
  80. py1(tim)=MSz(tim)./(((MSz(tim)).^2+(MSy(tim)).^2).^(0.5));
  81. pz1(tim)=-MSy(tim)./(((MSz(tim)).^2+(MSy(tim)).^2).^(0.5));
  82. px2(tim)=MSz(tim)./(((MSz(tim)).^2+(MSx(tim)).^2).^(0.5));
  83. py2(tim)=0;
  84. pz2(tim)=-MSx(tim)./(((MSz(tim)).^2+(MSx(tim)).^2).^(0.5));
  85. vx(tim)=(MSx(tim))./(MS(tim));
  86. vy(tim)=(MSy(tim))./(MS(tim));
  87. vz(tim)=(MSz(tim))./(MS(tim));
  88. cpx(tim)=MCXd(tim)+RM.*vx(tim)./(sin(MSa(tim)));
  89. cpy(tim)=MCYd(tim)+RM.*vy(tim)./(sin(MSa(tim)));
  90. cpz(tim)=MCZ(tim)+RM.*vz(tim)./(sin(MSa(tim)));
  91. opx(tim)=MCXd(tim)-RM.*vx(tim)./(sin(MSo(tim)));
  92. opy(tim)=MCYd(tim)-RM.*vy(tim)./(sin(MSo(tim)));
  93. opz(tim)=MCZ(tim)-RM.*vz(tim)./(sin(MSo(tim)));
  94. %vex='vex';vey='vey';vez='vez';%a='a';b='b';
  95. vex=zeros(length(ti),length(a),length(b));vey=vex;vez=vex;vox=vex;voy=vox;voz=vox;
  96. cosbMSa=repmat(cos(b.*MSa(index)),[length(a) 1]);
  97. cosbMSo=repmat(cos(b.*MSo(index)+(1-b).*MSa(index)),[length(a) 1]);
  98. try
  99. for index=1:length(ti)
  100.     vex(index,:,:)=(vx(index)).*cos(b.*MSa(index))+(px1(index)).*cos(a')*sin(b.*MSa(index))+(px2(index)).*sin(a')*sin(b.*MSa(index));
  101.     vey(index,:,:)=(vy(index)).*cos(b.*MSa(index))+(py1(index)).*cos(a')*sin(b.*MSa(index))+(py2(index)).*sin(a')*sin(b.*MSa(index));
  102.     vez(index,:,:)=(vz(index)).*cos(b.*MSa(index))+(pz1(index)).*cos(a')*sin(b.*MSa(index))+(pz2(index)).*sin(a')*sin(b.*MSa(index));
  103.     vox(index,:,:)=(vx(index)).*cos(b.*MSo(index)+(1-b).*MSa(index))+(px1(index)).*cos(a')*sin(b.*MSo(index)+(1-b).*MSa(index))+(px2(index)).*sin(a')*sin(b.*MSo(index)+(1-b).*MSa(index));
  104.     voy(index,:,:)=(vy(index)).*cos(b.*MSo(index)+(1-b).*MSa(index))+(py1(index)).*cos(a')*sin(b.*MSo(index)+(1-b).*MSa(index))+(py2(index)).*sin(a')*sin(b.*MSo(index)+(1-b).*MSa(index));
  105.     voz(index,:,:)=(vz(index)).*cos(b.*MSo(index)+(1-b).*MSa(index))+(pz1(index)).*cos(a')*sin(b.*MSo(index)+(1-b).*MSa(index))+(pz2(index)).*sin(a')*sin(b.*MSo(index)+(1-b).*MSa(index));
  106. end
  107. ck(tim,:,:)=(((((cpx(tim)').*vex(tim,:,:)+(cpy(tim)').*vey(tim,:,:))/(REa^2)+((cpz(tim)').*vez(tim,:,:))/(REb^2)).^2-(((vex(tim,:,:)).^2+(vey(tim,:,:)).^2)/(REa^2)+((vez(tim,:,:)).^2)/(REb^2)).*(((cpx(tim)').^2+(cpy(tim)').^2)/(REa.^2)+((cpz(tim)').^2)/(REb^2)-1)).^0.5)./(((vex(tim,:,:)).^2+(vey(tim,:,:)).^2)/(REa^2)+((vez(tim,:,:)).^2)/(REb^2))-(((cpx(tim)').*vex(tim,:,:)+(cpy(tim)').*vey(tim,:,:))/(REa.^2)+((cpz(tim)').*vez(tim,:,:))/(REb^2))./(((vex(tim,:,:)).^2+(vey(tim,:,:)).^2)/(REa^2)+((vez(tim,:,:)).^2)/(REb^2));
  108. cpxx(tim,:,:)=cpx(tim)'.*ones(length(ti),length(a),length(b))+ck(tim,:,:).*(vex(tim,:,:));
  109. cpyy(tim,:,:)=cpy(tim)'.*ones(length(ti),length(a),length(b))+ck(tim,:,:).*vey(tim,:,:);
  110. cpzz(tim,:,:)=cpz(tim)'.*ones(length(ti),length(a),length(b))+ck(tim,:,:).*vez(tim,:,:);
  111.  
  112. ok(tim,:,:)=((((((opx(tim)').*vox(tim,:,:)+(opy(tim)').*voy(tim,:,:))/(REa^2)+((opz(tim)').*voz(tim,:,:))/(REb^2)).^2-(((vox(tim,:,:)).^2+(voy(tim,:,:)).^2)/(REa^2)+((voz(tim,:,:)).^2)/(REb^2)).*(((opx(tim)').^2+(opy(tim)').^2)/(REa.^2)+((opz(tim)').^2)/(REb^2)-1)).^0.5)./(((vox(tim,:,:)).^2+(voy(tim,:,:)).^2)/(REa^2)+((voz(tim,:,:)).^2)/(REb^2))-(((opx(tim)').*vox(tim,:,:)+(opy(tim)').*voy(tim,:,:))/(REa.^2)+((opz(tim)').*voz(tim,:,:))/(REb^2))./(((vox(tim,:,:)).^2+(voy(tim,:,:)).^2)/(REa^2)+((voz(tim,:,:)).^2)/(REb^2)));
  113.  
  114. opxx(tim,:,:)=opx(tim)'.*ones(length(ti),length(a),length(b))+ok(tim,:,:).*(vox(tim,:,:));
  115. opyy(tim,:,:)=opy(tim)'.*ones(length(ti),length(a),length(b))+ok(tim,:,:).*voy(tim,:,:);
  116. opzz(tim,:,:)=opz(tim'.*ones(length(ti),length(a),length(b)))+ok(tim,:,:).*voz(tim,:,:);
  117. catch
  118. for index=1:length(ti)
  119.     vex(index,:,:)=(vx(index)).*cosbMSa+(px1(index)).*cos(a')*sin(b.*MSa(index))+(px2(index)).*sin(a')*sin(b.*MSa(index));
  120.     vey(index,:,:)=(vy(index)).*cosbMSa+(py1(index)).*cos(a')*sin(b.*MSa(index))+(py2(index)).*sin(a')*sin(b.*MSa(index));
  121.     vez(index,:,:)=(vz(index)).*cosbMSa+(pz1(index)).*cos(a')*sin(b.*MSa(index))+(pz2(index)).*sin(a')*sin(b.*MSa(index));
  122.     vox(index,:,:)=(vx(index)).*cosbMSo+(px1(index)).*cos(a')*sin(b.*MSo(index)+(1-b).*MSa(index))+(px2(index)).*sin(a')*sin(b.*MSo(index)+(1-b).*MSa(index));
  123.     voy(index,:,:)=(vy(index)).*cosbMSo+(py1(index)).*cos(a')*sin(b.*MSo(index)+(1-b).*MSa(index))+(py2(index)).*sin(a')*sin(b.*MSo(index)+(1-b).*MSa(index));
  124.     voz(index,:,:)=(vz(index)).*cosbMSo+(pz1(index)).*cos(a')*sin(b.*MSo(index)+(1-b).*MSa(index))+(pz2(index)).*sin(a')*sin(b.*MSo(index)+(1-b).*MSa(index));
  125. end
  126. cpxtim=repmat(cpx(tim)',[1 length(a) length(b)]);
  127. cpytim=repmat(cpy(tim)',[1 length(a) length(b)]);
  128. cpztim=repmat(cpz(tim)',[1 length(a) length(b)]);
  129. ck(tim,:,:)=((((cpxtim.*vex(tim,:,:)+cpytim.*vey(tim,:,:))/(REa^2)+(cpztim.*vez(tim,:,:))/(REb^2)).^2-(((vex(tim,:,:)).^2+(vey(tim,:,:)).^2)/(REa^2)+((vez(tim,:,:)).^2)/(REb^2)).*((cpxtim.^2+cpytim.^2)/(REa.^2)+(cpztim.^2)/(REb^2)-1)).^0.5)./(((vex(tim,:,:)).^2+(vey(tim,:,:)).^2)/(REa^2)+((vez(tim,:,:)).^2)/(REb^2))-((cpxtim.*vex(tim,:,:)+cpytim.*vey(tim,:,:))/(REa.^2)+(cpztim.*vez(tim,:,:))/(REb^2))./(((vex(tim,:,:)).^2+(vey(tim,:,:)).^2)/(REa^2)+((vez(tim,:,:)).^2)/(REb^2));
  130. cpxx(tim,:,:)=cpxtim+ck(tim,:,:).*(vex(tim,:,:));
  131. cpyy(tim,:,:)=cpytim+ck(tim,:,:).*vey(tim,:,:);
  132. cpzz(tim,:,:)=cpztim+ck(tim,:,:).*vez(tim,:,:);
  133.  
  134. opxtim=repmat(opx(tim)',[1 length(a) length(b)]);
  135. opytim=repmat(opy(tim)',[1 length(a) length(b)]);
  136. opztim=repmat(opz(tim)',[1 length(a) length(b)]);
  137.  
  138. ok(tim,:,:)=(((((opxtim.*vox(tim,:,:)+opytim.*voy(tim,:,:))/(REa^2)+(opztim.*voz(tim,:,:))/(REb^2)).^2-(((vox(tim,:,:)).^2+(voy(tim,:,:)).^2)/(REa^2)+((voz(tim,:,:)).^2)/(REb^2)).*((opxtim.^2+opytim.^2)/(REa.^2)+(opztim.^2)/(REb^2)-1)).^0.5)./(((vox(tim,:,:)).^2+(voy(tim,:,:)).^2)/(REa^2)+((voz(tim,:,:)).^2)/(REb^2))-((opxtim.*vox(tim,:,:)+opytim.*voy(tim,:,:))/(REa.^2)+(opztim.*voz(tim,:,:))/(REb^2))./(((vox(tim,:,:)).^2+(voy(tim,:,:)).^2)/(REa^2)+((voz(tim,:,:)).^2)/(REb^2)));
  139.  
  140. opxx(tim,:,:)=opxtim+ok(tim,:,:).*(vox(tim,:,:));
  141. opyy(tim,:,:)=opytim+ok(tim,:,:).*voy(tim,:,:);
  142. opzz(tim,:,:)=opztim+ok(tim,:,:).*voz(tim,:,:);
  143. end
  144.  
  145.  
  146. %Sun shadow - night time - on Earth
  147. zxya(tim)=atan2((SCZ(tim)),(((SCXd(tim)).^2+(SCYd(tim)).^2).^0.5));
  148. xya(tim)=atan2((SCYd(tim)),(SCXd(tim)));
  149. adda(tim,1:length(b))=asin((tan(zxya(tim)'))*(tan(b*pi+0.5*pi))*(REa/REb));
  150. size(adda)
  151.  
  152. %Maple just discards complex parts of numbers, unlike Matlab
  153. %REa*cos(b*Pi+0.5*Pi)*cos(0.5*a+((adda(ti))*(1-(a/Pi)))+0.5*Pi+xya(ti)),REa*cos(b*Pi+0.5*Pi)*sin(0.5*a+((adda(ti))*(1-a/Pi))+0.5*Pi+xya(ti)),REb*sin(b*Pi+0.5*Pi)
  154. %plot3(REa*cos(b*pi+0.5*pi)*cos(0.5*ones(length(ti),1)*a+(reshape(adda,[2,length(ti),length(b),1])*(1-(a/pi)))+0.5*pi+xya(:,tim)),REa.*cos(b.*pi+0.5.*pi).*sin(0.5.*a+((adda(tim)).*(1-a/pi))+0.5.*pi+xya(tim)),REb.*sin(b.*pi+0.5.*pi),'m');
  155.  
  156. try
  157.     colormap winter;
  158. catch
  159.    
  160. end
  161. figure;
  162. s=surf(REa*cos(a')*cos(b*2*pi),REa*cos(a')*sin(b*2*pi),REb*sin(a')*ones(1,length(b)));%Earth itself
  163. try
  164.     s.EdgeColor = 'none';
  165. catch
  166.     set(s,'color','b')
  167. end
  168.  
  169. hold on;%If done before surf, hold on makes figure 2D
  170.  
  171.  xlim([-7 7]) %Focusing on Earth
  172.  ylim([-7 7])
  173.  zlim([-7 7])
  174.  
  175. for andex=1:length(a)
  176.     addapi(tim,1:length(b),andex)=((1-(a(andex)/pi))*(adda(tim,:)))+0.5*a(andex)+0.5*pi;
  177. end
  178. try
  179.     addapi(tim,:,:,:)=addapi(tim,:,:,:)+xya(:,tim)';
  180. catch
  181.     addapi(tim,:,:,:)=addapi(tim,:,:,:)+repmat(xya(:,tim)',[1 length(b) length(a)]);
  182. end
  183. for andex=1:length(a)
  184.     for bndex=1:length(b)
  185.         darkx(tim,andex,bndex)=REa*cos(b(bndex)*pi+0.5*pi)*cos(addapi(tim,bndex,andex));
  186.         darky(tim,andex,bndex)=REa*cos(b(bndex)*pi+0.5*pi)*sin(addapi(tim,bndex,andex));
  187.         darkz(tim,andex,bndex)=REb*sin(b(bndex)*pi+0.5*pi);
  188.     end
  189. end
  190. start=round(0.5*length(ti));%1
  191. night=surf(squeeze(real(darkx(start,:,:))),squeeze(real(darky(start,:,:))),squeeze(real(darkz(start,:,:))));
  192. try
  193.     night.EdgeColor = 'black';
  194. catch
  195. end
  196. %hold on;
  197. c=surf(squeeze(real(cpxx(start,:,:))),squeeze(real(cpyy(start,:,:))),squeeze(real(cpzz(start,:,:))));
  198. try
  199.     c.EdgeColor = 'yellow';
  200. catch
  201. end
  202. %hold on;
  203. o=surf(squeeze(real(opxx(start,:,:))),squeeze(real(opyy(start,:,:))),squeeze(real(opzz(start,:,:))));
  204. try
  205.     o.EdgeColor = 'red';
  206. catch
  207. end
  208. %hold on;
  209. try
  210. for index=1:length(ti)
  211.     set(night,'XData',squeeze(real(darkx(index,:,:))));
  212.     set(night,'YData',squeeze(real(darky(index,:,:))));
  213.     set(c,'XData',squeeze(real(cpxx(index,:,:))));
  214.     set(c,'YData',squeeze(real(cpyy(index,:,:))));
  215.     set(c,'ZData',squeeze(real(cpzz(index,:,:))));
  216.     set(o,'XData',squeeze(real(opxx(index,:,:))));
  217.     set(o,'YData',squeeze(real(opyy(index,:,:))));
  218.     set(o,'ZData',squeeze(real(opzz(index,:,:))));
  219.     pause(1);
  220. end
  221. catch
  222.     Ex=REa*cos(a')*cos(b*2*pi);
  223.     Ey=REa*cos(a')*sin(b*2*pi);
  224.     Ez=REb*sin(a')*ones(1,length(b));
  225.     for index=min(13,length(ti)):max(13,length(ti)-13)
  226.         clf;
  227.         start=index;
  228.         nightx=squeeze(real(darkx(start,:,:)));
  229.         nighty=squeeze(real(darky(start,:,:)));
  230.         nightz=squeeze(real(darkz(start,:,:)));
  231.         cx=squeeze(real(cpxx(start,:,:)));
  232.         cy=squeeze(real(cpyy(start,:,:)));
  233.         cz=squeeze(real(cpzz(start,:,:)));
  234.         ox=squeeze(real(opxx(start,:,:)));
  235.         oy=squeeze(real(opyy(start,:,:)));
  236.         oz=squeeze(real(opzz(start,:,:)));
  237.     %s=surf([Ex,Ey,Ez],[nightx,nighty,nightz],[cx,cy,cz],[ox,oy,oz]);%everything
  238.         s=surf([Ex,nightx,cx,ox],[Ey,nighty,cy,oy],[Ez,nightz,cz,oz]);%everything
  239.         pause(0.2);
  240.     end
  241. end
  242. %plot3(REa*cos(a')*cos(b*2*pi),REa*cos(a')*sin(b*2*pi),REb*sin(a')*ones(1,length(b)),'b');
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!