[IDL] 克卜勒行星運動定律、相差日距程式

這是一份尚未完成的程式，不過我暫時要封存它十天(直到我下週三 26號報告完…)
所以先把它上傳這邊放著…



這是傅老的太陽系(大二選修課)的作業
身為助教(汗)，應該要能夠從最基本的層面理解這個計算的內容…
不過我真的是太弱了… 要從高中數學開始慢慢理解

中間幾經波折，最大的問題在得到五個火星位置座標後，如何得出橢圓軌道方程？
最終的突破就在一小時前，我們松山的物理實習老師翁韶君(她坐在我後面)帶來的書完美的提供了我所需的一切：
『物理發展史』授課教授：姚珩(師大物理系)

姚教授用天文的角度講解物理發展史，著重(幾乎整門課)於克卜勒和牛頓定律的推導
這種紮實的演算和簡單概念的理解，是我現在看來非常、非常重要的
以前在學習天文時，總覺得古人得出的結論和公式理所當然
雖然知道背後的過程肯定艱辛，但公式是優雅的，就不覺得有什麼需要證明的
然後不知不覺中，就在學習的物理和數學(相對基礎紮實)和天文學間產生了隱形的鴻溝
如何從簡單的橢圓方程和觀測資料得出行星運動定律？
尤其，克卜勒根本沒有使用微積分(還沒發明)，等同於高中程度的我們，應該也能做出同樣的發現才是

對我而且這是非常薄弱的… 我很多學生學弟妹都比我強很多，我要努力才能補齊
也難怪傅老總是對我們那麼弱感到生氣… 老師我對不起你!!

以下是程式碼，反正總是會更新，是個隨便寫寫，沒頭沒尾的程式XD
(啊對了，這份程式是我第一個實際應用function功能的程式… 早在太陽model就該會用了T^T)






;This program is for 2011/10/04 Class:SolorSystem's Homework.
;Data from HHF with date, Earth heliocentric longitude, Mars geocentric longitude
;Written by Hsieh Hsiang-yu. Update:2011/10/17

function year2date1583A,yr,mt,dt
  if yr ge 1592 then dt=dt+1.
    if yr eq 1592 and mt le 2 then dt=dt-1.
  if yr ge 1588 then dt=dt+1.
    if yr eq 1588 and mt le 2 then dt=dt-1.
  if yr ge 1584 then dt=dt+1.
    if yr eq 1584 and mt le 2 then dt=dt-1.
  if yr gt 1583 then dt=dt+285.
bmt=31.
smt=30.
Fmt=28.
  if yr eq 1583 and mt eq 3 and dt lt 21 then dt=double(dt-21.)
  if yr eq 1583 and mt eq 2 then dt=double(dt-21.-fmt)
  if yr eq 1583 and mt eq 1 then dt=double(dt-21.-fmt-bmt)
  if yr eq 1583 then yr=yr+1.
  if mt eq 1 then mt=double(0.)
  if mt eq 2 then mt=double(bmt)
  if mt eq 3 then mt=double(bmt+Fmt)
  if mt eq 4 then mt=double(2*bmt+Fmt)
  if mt eq 5 then mt=double(2*bmt+Fmt+smt)
  if mt eq 6 then mt=double(3*bmt+Fmt+smt)
  if mt eq 7 then mt=double(3*bmt+Fmt+2*smt)
  if mt eq 8 then mt=double(4*bmt+Fmt+2*smt)
  if mt eq 9 then mt=double(5*bmt+Fmt+2*smt)
  if mt eq 10 then mt=double(5*bmt+Fmt+3*smt)
  if mt eq 11 then mt=double(6*bmt+Fmt+3*smt)
  if mt eq 12 then mt=double(6*bmt+Fmt+4*smt)
Return,double((yr-1-1583)*365+mt+dt)
end

pro KeplerMarsOrbit
close,/all
print,'====start!===='
Q='  '
StartDate=double(1583.+(31.+28.+21.)/365.)
openr,1,'D:\Kepler\MarsData.txt'
readf,1,Q
dN=10
dR=dblarr(dN)
EH=dblarr(dN)
MG=dblarr(dN)
EPX=dblarr(dN)
EPY=dblarr(dN)

for i=0,9 do begin
readf,1,format='(f4,2(4x,f2),2(4x,f3,x,f2))',yr,mt,dt,EHd,EHm,MGd,MGm
;dR(i)=year2date1583A(yr,mt,dt)
;if dR(i) gt 365.2422 then dR(i)=dR(i)-fix(dR(i)/365.2422)*365.2422
EH(i)=double(EHd+EHm/60.0d0)
MG(i)=double(MGd+MGm/60.0d0)
EPX(i)=cos(EH(i)*!dtor)
EPY(i)=sin(EH(i)*!dtor)
endfor

set_plot,'ps'
device, filename="D:\Kepler\MarsOrbit.ps",/color
TVLCT,[0,255,0,0,255],[0,0,255,0,255],[0,0,0,255,0] ;0=black,1=red,2=green,3=blue,4=yellow
A=FINDGEN(33) * (!PI*2/32.)
USERSYM, COS(A), SIN(A);, /FILL
plot,EPX,EPY,psym=8,xtitle='(AU)',ytitle='(AU)',title='Earth & Mars Orbit',$
      symsize=1,xrange=[-2,2],yrange=[-2,2],/isotropic,/xstyle,/ystyle,color=0
oplot,COS(A),SIN(A),psym=0,color=0
oplot,0*EPX,0*EPY,psym=1

for i=0,9 do begin
B=findgen(10)/4.
E2MX=EPX(i)+B*(cos(MG(i)*!dtor))
E2MY=EPY(i)+B*(sin(MG(i)*!dtor))
oplot,E2MX,E2MY,psym=0,color=1
endfor

MPX=dblarr(dN/2)
MPY=dblarr(dN/2)
u=0
for i=0,4 do begin
ii=(sin(MG(i*2+1)*!dtor)*(EPX(i*2+1)-EPX(i*2))-cos(MG(i*2+1)*!dtor)*(EPY(i*2+1)-EPY(i*2)))/$
       (cos(MG(i*2)*!dtor)*sin(MG(i*2+1)*!dtor)-sin(MG(i*2)*!dtor)*cos(MG(i*2+1)*!dtor))
;jj=(sin(MG(i*2)*!dtor)*(EPX(i*2)-EPX(i*2+1))-cos(MG(i*2)*!dtor)*(EPY(i*2)-EPY(i*2+1)))/$
;       (cos(MG(i*2+1)*!dtor)*sin(MG(i*2)*!dtor)-sin(MG(i*2+1)*!dtor)*cos(MG(i*2)*!dtor))
MPX(i)=EPX(i*2)+cos(MG(i*2)*!dtor)*(ii)
MPY(i)=EPY(i*2)+sin(MG(i*2)*!dtor)*(ii)
endfor
oplot,MPX,MPY,psym=8,color=1

RMS=double(SQRT(MPX^2.+MPY^2.))
RMS=double(1./RMS)
DMS=dblarr(5)
for i=0,4 do begin
  DMS(i)=atan(MPY(i)/MPX(i))
  if MPX(i) lt 0. and MPY(i) gt 0. then DMS(i)=DMS(i)-!dpi  ;2-quarter
  if MPX(i) lt 0. and MPY(i) lt 0. then DMS(i)=DMS(i)-!dpi  ;3-quarter
  if MPX(i) gt 0. and MPY(i) lt 0. then DMS(i)=DMS(i)+2*!dpi;4-quarter
endfor
for i=0,4 do begin
  B=findgen(11)/(RMS(i)*10.)
  S2MX=B*cos(DMS(i))
  S2MY=B*sin(DMS(i))
  oplot,S2MX,S2MY,psym=0,color=2 ;0=black,1=red,2=green,3=blue,4=yellow
endfor

DDMS=[transpose(cos(DMS)),transpose(sin(DMS))]
result=regress(DDMS,RMS,SIGMA=sss,CONST=ccc,MEASURE_ERRORS=mmm,/double)
  ee=(sqrt(result[0]^2+result[1]^2))/ccc
  print,ee

device,/close
close,/all
print,'====DONE!===='
end

=====================================

我要快點把論文要做的東西弄好… 但還是一句老話，要看書啊…T^T