Outline

plan de l’ecliptique

completely misleading

path of sun across the sky

the sun path along the ecliptic through the zodiac

date ?S ?M r 15/05/1992,55,21,7,44,0.7324 2/04/1994,12?56’,350?48’,0,7241 18/02/1996,329?29’,332?41’,0,7156 5/01/1998,285?16’,314?34’,0,7122 23/11/1999,241?01’,298?05’,0,7150 10/10/2001,197?27’,288?40’,0,7231 15/07/2005,113?17’,22?22’,0,7362

getwd()
## [1] "C:/Users/phili/OneDrive/Desktop/ForBOOK3450/Pres-KeplerMars"

Reconstruction of the orbit of the Earth

mydata=read.csv("keplermars1.csv")
#summary(mydata)
head(mydata)
##         date lsdeg lsmin lmdeg lmmin      r
## 1 15/05/1992    55    21     7    44 0.7324
## 2  2/04/1994    12    56   350    48 0.7241
## 3 18/02/1996   329    29   332    41 0.7156
## 4  5/01/1998   285    16   314    34 0.7122
## 5 23/11/1999   241     1   298     5 0.7150
## 6 10/10/2001   197    27   288    40 0.7231
ls = pi *(mydata$lsdeg + mydata$lsmin/60)/180.
lm = pi *(mydata$lmdeg + mydata$lmmin/60)/180.

ls
## [1] 0.9660397 0.2257292 5.7505690 4.9788426 4.2065344 3.4461526 1.9771672
lm
## [1] 0.1349721 6.1226150 5.8064195 5.4902241 5.2025356 5.0381838 0.3903720
lmdeg0=334  
lmmin0=58
lsdeg0 =154 
lsmin0 =58

ls0=pi *(lsdeg0 + lsmin0/60)/180.
lm0=pi *(lmdeg0 + lmmin0/60)/180.
xd=c()
yd=c()
SM0=1
for(i in 1:7){
theta=ls[i]-lm0+pi
r = SM0 *  sin(lm[i]-lm0) / sin(ls[i]-lm[i]) 
xd=c(xd,r * cos(theta))
yd=c(yd,r * sin(theta))
}
xd
## [1] -0.12240368 -0.57091171 -0.71045783 -0.46068844  0.04926044  0.53328958
## [7]  0.54981097
yd
## [1] -0.72241597 -0.44551047  0.06820073  0.54322516  0.71340309  0.48838456
## [7] -0.48957689
XY <- matrix(c(xd,yd),7,2,byrow=FALSE)
XY
##             [,1]        [,2]
## [1,] -0.12240368 -0.72241597
## [2,] -0.57091171 -0.44551047
## [3,] -0.71045783  0.06820073
## [4,] -0.46068844  0.54322516
## [5,]  0.04926044  0.71340309
## [6,]  0.53328958  0.48838456
## [7,]  0.54981097 -0.48957689
plot(XY,xlim=c(-1,1),ylim=c(-1,1),asp=1)

sqrt(xd^2+yd^2)
## [1] 0.7327124 0.7241683 0.7137238 0.7122692 0.7151018 0.7231302 0.7361913
library(conicfit)
## Warning: package 'conicfit' was built under R version 3.6.3
## Loading required package: pracma
## Loading required package: geigen
xd
## [1] -0.12240368 -0.57091171 -0.71045783 -0.46068844  0.04926044  0.53328958
## [7]  0.54981097
yd
## [1] -0.72241597 -0.44551047  0.06820073  0.54322516  0.71340309  0.48838456
## [7] -0.48957689
xy<-calculateEllipse(0,0,200,100,45,50, randomDist=TRUE,noiseFun=function(x)
(x+rnorm(1,mean=0,sd=50)))
plot(xy[,1],xy[,2],xlim=c(-250,250),ylim=c(-250,250),col='magenta');par(new=TRUE)
ellipDirect <- EllipseDirectFit(xy)
ellipDirectG <- AtoG(ellipDirect)$ParG
xyDirect<-calculateEllipse(ellipDirectG[1], ellipDirectG[2], ellipDirectG[3],
ellipDirectG[4], 180/pi*ellipDirectG[5])
plot(xyDirect[,1],xyDirect[,2],xlim=c(-250,250),ylim=c(-250,250),type='l',
col='cyan');par(new=TRUE)

xy
##               [,1]        [,2]
##  [1,] -184.6270161  -50.138878
##  [2,]   16.7471001   28.343658
##  [3,] -163.1755614   -5.868783
##  [4,]  -50.6315612 -180.117007
##  [5,]   97.5354061   33.566471
##  [6,] -181.5608039 -110.320917
##  [7,]   10.2394705   23.856691
##  [8,]  -90.0217623   87.712313
##  [9,]  -72.1310219 -187.261003
## [10,]  185.0545180   88.761389
## [11,]    3.3164010  -95.638459
## [12,]  156.4003131  134.004194
## [13,]  216.0506752   12.553730
## [14,]  -19.0633433  133.836828
## [15,]  110.6943296  135.741257
## [16,]  -66.8306982   76.839398
## [17,]  121.4220890   79.444007
## [18,]   66.1150237  141.013513
## [19,] -139.5438726  -26.034570
## [20,] -170.6545342 -160.855429
## [21,] -103.9736774 -165.570715
## [22,]  -20.6123674 -140.058205
## [23,]  152.4195871  -71.857464
## [24,] -150.8087287 -139.482252
## [25,]  -37.8305628 -140.332590
## [26,] -150.3226151  -87.958621
## [27,]  151.2711175   49.755999
## [28,] -132.9095328 -126.989811
## [29,] -150.4949728 -186.409207
## [30,]   92.9446615  -76.815457
## [31,] -140.9885746   -2.683170
## [32,]  155.4394629   46.348636
## [33,]    6.5913582 -174.004220
## [34,]  121.7970490   94.203069
## [35,]   65.5364053  141.715146
## [36,]  -37.3055671  -67.038866
## [37,]   51.8938538  237.496855
## [38,] -151.9446468  -95.946510
## [39,]   50.2423218  -32.284310
## [40,]  -80.4634526 -177.862618
## [41,]   -0.2358382  141.560851
## [42,] -194.6634251 -118.039037
## [43,]   65.3909840  239.825208
## [44,]   44.3938268 -171.410242
## [45,]   69.5293981    8.187235
## [46,] -170.2222705 -147.139124
## [47,]  -45.8206400   -9.667809
## [48,]  183.0922605  115.623779
## [49,]  -64.3357786 -219.324133
## [50,] -101.8555933 -183.135565
XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=TRUE)
XY
##      [,1] [,2]
## [1,]    1    7
## [2,]    2    6
## [3,]    5    8
## [4,]    7    7
## [5,]    9    5
## [6,]    3    7
## [7,]    6    2
## [8,]    8    4
XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=FALSE)
XY
##      [,1] [,2]
## [1,]    1    9
## [2,]    7    5
## [3,]    2    3
## [4,]    6    7
## [5,]    5    6
## [6,]    8    2
## [7,]    7    8
## [8,]    7    4
XY <- matrix(c(xd,yd),7,2,byrow=FALSE)
XY
##             [,1]        [,2]
## [1,] -0.12240368 -0.72241597
## [2,] -0.57091171 -0.44551047
## [3,] -0.71045783  0.06820073
## [4,] -0.46068844  0.54322516
## [5,]  0.04926044  0.71340309
## [6,]  0.53328958  0.48838456
## [7,]  0.54981097 -0.48957689
plot(XY)

sqrt(xd^2+yd^2)
## [1] 0.7327124 0.7241683 0.7137238 0.7122692 0.7151018 0.7231302 0.7361913
library(conicfit)
#xy=c(xd,yd)

XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=TRUE)
XY
##      [,1] [,2]
## [1,]    1    7
## [2,]    2    6
## [3,]    5    8
## [4,]    7    7
## [5,]    9    5
## [6,]    3    7
## [7,]    6    2
## [8,]    8    4
XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=FALSE)
XY
##      [,1] [,2]
## [1,]    1    9
## [2,]    7    5
## [3,]    2    3
## [4,]    6    7
## [5,]    5    6
## [6,]    8    2
## [7,]    7    8
## [8,]    7    4
XY <- matrix(c(xd,yd),7,2,byrow=FALSE)
XY
##             [,1]        [,2]
## [1,] -0.12240368 -0.72241597
## [2,] -0.57091171 -0.44551047
## [3,] -0.71045783  0.06820073
## [4,] -0.46068844  0.54322516
## [5,]  0.04926044  0.71340309
## [6,]  0.53328958  0.48838456
## [7,]  0.54981097 -0.48957689
plot(XY)

xy=XY


#EllipseDirectFit(XY)

ellipDirect <- EllipseDirectFit(xy)
ellipDirectG <- AtoG(ellipDirect)$ParG
xyDirect<-calculateEllipse(ellipDirectG[1], ellipDirectG[2], ellipDirectG[3],
ellipDirectG[4], 180/pi*ellipDirectG[5])
plot(xyDirect[,1],xyDirect[,2],type='l',
col='cyan',asp=1);par(new=TRUE)
lines(XY)

ellipDirect
##               [,1]
## [1,]  0.6644022619
## [2,] -0.0001387037
## [3,]  0.6615275858
## [4,] -0.0103261978
## [5,]  0.0129203166
## [6,] -0.3473863892
print('a')
## [1] "a"
b=-ellipDirect[2]/ellipDirect[1]
a=-ellipDirect[3]/ellipDirect[1]
c=-ellipDirect[4]/ellipDirect[1]
d=-ellipDirect[5]/ellipDirect[1]
f=-ellipDirect[6]/ellipDirect[1]

a
## [1] -0.9956733
b
## [1] 0.0002087646
c
## [1] 0.01554209
d
## [1] -0.01944653
f
## [1] 0.5228555
-xd^2+a*yd^2+b*xd*yd+c*xd+d*yd+f
## [1]  0.0004106042 -0.0008619026  0.0010955928 -0.0009712575  0.0005867508
## [6] -0.0001843060 -0.0000754816
gett1 <-function(i,s){
print(c(i,s,ls1[i],lm1[i]))
  
theta1 = ls1[i]-lm0+pi
t=tan(theta1)
aa=a*t^2 + b*t-1
bb=c+d*t
cc=f
delta=bb*bb-4*aa*cc

x1=(-bb+s*sqrt(delta))/(2*aa)
y1=x1*t
return(c(x1,y1))
}

gett2 <-function(i,s){
print(c(i,s,ls2[i],lm2[i]))
  
theta2 = ls2[i]-lm0+pi
t=tan(theta2)
aa=a*t^2 + b*t-1
bb=c+d*t
cc=f
delta=bb*bb-4*aa*cc

x2=(-bb+s*sqrt(delta))/(2*aa)
y2=x2*t
return(c(x2,y2))
}

[1,] 0.6644022619 [2,] -0.0001387037 [3,] 0.6615275858 [4,] -0.0103261978 [5,] 0.0129203166 [6,] -0.3473863892

“a” [1] -0.9956733 [1] 0.0002087646 [1] 0.01554209 [1] -0.01944653 [1] 0.5228555

Beatrice Sandre

a=-0.9996 b=-0.0002594 c=0.01495 d=-0.01926 f=0.5242

XY <- matrix(c(xd,yd),7,2,byrow=FALSE)
xy=XY

ellipTaubin <- EllipseFitByTaubin(xy)
ellipTaubinG <- AtoG(ellipTaubin)$ParG
xyTaubin<-calculateEllipse(ellipTaubinG[1], ellipTaubinG[2], ellipTaubinG[3],
ellipTaubinG[4], 180/pi*ellipTaubinG[5])
plot(xyTaubin[,1],xyTaubin[,2],xlim=c(-1,1),ylim=c(-1,1),type='l',
col='red');par(new=TRUE)

lines(XY)

ellipTaubin
##               [,1]
## [1,] -0.6644018823
## [2,]  0.0001386031
## [3,] -0.6615279421
## [4,]  0.0103264335
## [5,] -0.0129203515
## [6,]  0.3473864285
print('a')
## [1] "a"
b=-ellipTaubin[2]/ellipTaubin[1]
a=-ellipTaubin[3]/ellipTaubin[1]
c=-ellipTaubin[4]/ellipTaubin[1]
d=-ellipTaubin[5]/ellipTaubin[1]
f=-ellipTaubin[6]/ellipTaubin[1]

a
## [1] -0.9956744
b
## [1] 0.0002086134
c
## [1] 0.01554245
d
## [1] -0.01944659
f
## [1] 0.5228559
xc=(2*a*c-b*d)/(4*a+b*b)
xc
## [1] 0.007770207
yc=(b*c+2*d)/(-b*b-4*a)
yc
## [1] -0.009764723
eps=f+xc*xc-a*yc*yc-b*xc*yc
eps
## [1] 0.5230112
phi=atan(yc/xc)

# longitude of perigee

lonp=180.*(phi+lm0)/pi
lonp     # 283 IMCCE
## [1] 283.4774
# coord of perigee

aa=a*(yc/xc)^2+b*(yc/xc)-1
bb=c+d*(yc/xc)
cc=f


print('delta')
## [1] "delta"
delta=bb-4*aa*cc
delta
## [1] 5.420575
xp=(-bb+sqrt(delta))/(2*aa)
yp=(yc/xc)*xp

# eccentricity
d1=sqrt(xc*xc+yc*yc)
d2=sqrt((xp-xc)^2+(yp-yc)^2)
e=d1/d2
e # 0.0168
## [1] 0.01717229
yc/xc
## [1] -1.256688
A=d2 # CP  demi grand axe 0.72
Ainv=1/d2 # IMCC 1.38

Reconstruction of the orbit of the Mars

mydata2=read.csv("keplermars2.csv")
#summary(mydata2)
mydata2
##        date1 ls1deg ls1min lm1deg lm1min      date2 ls2deg ls2min lm2deg lm2min
## 1 01/01/2003    280     54    230      2 18/11/2004    236     40    214     57
## 2 01/04/2003     11     36    287     34 16/02/2005    328      7    277      1
## 3 01/07/2003     99     29    335     26 18/05/2005     57     51    342     42
## 4 01/10/2003    188      8    330     12 18/08/2005    145     51     41     25
## 5 01/01/2002    281     10    347     26 19/11/2003    236     58    345     29
## 6 01/04/2002     11     50     51     47 17/02/2004    328     21     39      5
## 7 01/07/2002     99     44    112     20 18/05/2004     58      3     97      9
## 8 01/10/2002    188     23    171      5 18/08/2004    146      3    155     13
ls1 = pi *(mydata2$ls1deg + mydata2$ls1min/60)/180.
ls2 = pi *(mydata2$ls2deg + mydata2$ls2min/60)/180.
lm1 = pi *(mydata2$lm1deg + mydata2$lm1min/60)/180.
lm2 = pi *(mydata2$lm2deg + mydata2$lm2min/60)/180.

date1,ls1deg,ls1min,lm1deg,lm1min,date2,ls1deg,ls1min,lm1deg,lm1min 01/01/2003,280,54,230,2,18/11/2004,236,40,214,57 01/04/2003,11,36,287,34,16/02/2005,328,7,277,1 01/07/2003,99,29,335,26,18/05/2005,57,51,342,42 01/10/2003,188,8,330,12,18/08/2005,145,51,41,25 01/01/2002,281,10,347,26,19/11/2003,236,58,345,29 01/04/2002,11,50,51,47,17/02/2004,328,21,39,5 01/07/2002,99,44,112,20,18/05/2004,58,3,97,9 01/10/2002,188,23,171,05,18/08/2004,146,3,155,13

Given two dates (January 1 2003 and November 18 2004), we have one matching position for Mars (\((x,y)\)) seen from two positions of the Earth (\((x_1,y_1),(x_2,y_2)\)).

Geometry tells us that:

\[y=y_1+(x-x_1)\tan(\lambda_{M_1} - \lambda_{M_0} )\] \[y=y_2+(x-x_2)\tan(\lambda_{M_2} - \lambda_{M_0} )\]

xm=c()
ym=c()
ii=1
if(ii == 1){

ls1deg=280
ls1min=54
ls2deg=236
ls2min=40
  
lm1deg=230
lm1min=2 
lm2deg=214
lm2min=57
}


if(ii == 2){
ls1deg=11
ls1min=36
ls2deg=328
ls2min=7
  
lm1deg=287
lm1min=34 
lm2deg=277
lm2min=1
}

if(ii == 3){
ls1deg=99
ls1min=29
ls2deg=57
ls2min=51
  
lm1deg=335
lm1min=26 
lm2deg=342
lm2min=42
}

if(ii == 4){
ls1deg=99
ls1min=29
ls2deg=57
ls2min=51
  
lm1deg=335
lm1min=26 
lm2deg=342
lm2min=42
}

if(ii == 5){
ls1deg=99
ls1min=29
ls2deg=57
ls2min=51
  
lm1deg=335
lm1min=26 
lm2deg=342
lm2min=42
}

if(ii == 6){
ls1deg=99
ls1min=29
ls2deg=57
ls2min=51
  
lm1deg=335
lm1min=26 
lm2deg=342
lm2min=42
}

if(ii == 7){
ls1deg=99
ls1min=29
ls2deg=57
ls2min=51
  
lm1deg=335
lm1min=26 
lm2deg=342
lm2min=42
}

if(ii == 8){
ls1deg=99
ls1min=29
ls2deg=57
ls2min=51
  
lm1deg=335
lm1min=26 
lm2deg=342
lm2min=42
}









ls1 = pi *(ls1deg + ls1min/60)/180.
ls2 = pi *(ls2deg + ls2min/60)/180.
lm1 = pi *(lm1deg + lm1min/60)/180.
lm2 = pi *(lm2deg + lm2min/60)/180.

To obtain the coordinates of Earth at a date

getxy <-function(){
  
dmi=10000
theta1 = ls1-lm0+pi
aa=a*(tan(theta1))^2 + b*tan(theta1)-1
bb=c+d*tan(theta1)
cc=f
print('delta1')
delta=bb-4*aa*cc
delta
x1a=(-bb+sqrt(delta))/(2*aa)
y1a=x1a*tan(theta1)
print('x1-a')
#x1a
#y1a


dse=sqrt(x1a^2+y1a^2)
print('x1a: dse')
print(dse)
print(-x1a^2+a*y1a^2+b*x1a*y1a+c*x1a+d*y1a+f)

dto=abs(0.72-dse)
if(dto <= dmi){
  x1=x1a
  y1=y1a
  dmi=dto
}

x1b=(-bb-sqrt(delta))/(2*aa)
y1b=x1b*tan(theta1)
print('x1-b')
x1b
y1b
dse=sqrt(x1b^2+y1b^2)
print('x1b: dse')
print(dse)
print(-x1b^2+a*y1b^2+b*x1b*y1b+c*x1b+d*y1b+f)
dto=abs(0.72-dse)
if(dto <= dmi){
  x1=x1b
  y1=y1b
  dmi=dto
}

dmi=10000
theta2 = ls2-lm0+pi
aa=a*(tan(theta2))^2 + b*tan(theta2)-1
bb=c+d*tan(theta2)
cc=f
print('delta2')
delta=bb-4*aa*cc
delta
print('x2-a')
x2a=(-bb+sqrt(delta))/(2*aa)
y2a=x2a*tan(theta2)

dse=sqrt(x2a^2+y2a^2)
print('x2a: dse')
print(dse)
print(-x2a^2+a*y2a^2+b*x2a*y2a+c*x2a+d*y2a+f)
dto=abs(0.72-dse)
if(dto <= dmi){
  x2=x2a
  y2=y2a
  dmi=dto
}

x2a
y2a
sqrt(x2a^2+y2a^2)
print('x2-b')
x2b=(-bb-sqrt(delta))/(2*aa)
y2b=x2b*tan(theta2)
dse=sqrt(x2b^2+y2b^2)
print('x2b: dse')
print(dse)
print(-x2b^2+a*y2b^2+b*x2b*y2b+c*x2b+d*y2b+f)
dto=abs(0.72-dse)
if(dto <= dmi){
  x2=x2b
  y2=y2b
  dmi=dto
}
x2b
y2b
sqrt(x2b^2+y2b^2)

print('rayon1 ST1')
print(sqrt(x1^2+y1^2))
print('rayon2 ST2')
print(sqrt(x2^2+y2^2))

#ax x + ay y = az
#bx x + by y = bz 

ax=tan(lm1-lm0)
ay=-1
az=-y1+x1*tan(lm1-lm0)
bx=tan(lm2-lm0)
by=-1
bz=-y2+x2*tan(lm2-lm0)
m=matrix(c(ax,ay,bx,by),nrow=2,byrow=TRUE)
mi=solve(m)
rhs=matrix(c(az,bz),nrow=2,byrow=TRUE)
vsol=mi%*%rhs

c(ax*vsol[1]+ay*vsol[2],az)
c(bx*vsol[1]+by*vsol[2],bz)

print(vsol)

vsoln=sqrt(vsol[1]^2+vsol[2]^2)

print(vsoln)

return(vsol)
}

vv=getxy()
## [1] "delta1"
## [1] "x1-a"
## [1] "x1a: dse"
## [1] 0.7141403
## [1] -0.003503242
## [1] "x1-b"
## [1] "x1b: dse"
## [1] 0.7390755
## [1] -0.003503242
## [1] "delta2"
## [1] "x2-a"
## [1] "x2a: dse"
## [1] 0.7327437
## [1] 0.0006886679
## [1] "x2-b"
## [1] "x2b: dse"
## [1] 0.7156717
## [1] 0.0006886679
## [1] "rayon1 ST1"
## [1] 0.7141403
## [1] "rayon2 ST2"
## [1] 0.7156717
##            [,1]
## [1,] -0.8026600
## [2,] -0.8599597
## [1] 1.176348
xm=c(xm,vv[1])
ym=c(ym,vv[2])
xm
## [1] -0.80266
ym
## [1] -0.8599597
xm
## [1] -0.80266
ym
## [1] -0.8599597
mydata2=read.csv("keplermars2.csv")
#summary(mydata2)
mydata2
##        date1 ls1deg ls1min lm1deg lm1min      date2 ls2deg ls2min lm2deg lm2min
## 1 01/01/2003    280     54    230      2 18/11/2004    236     40    214     57
## 2 01/04/2003     11     36    287     34 16/02/2005    328      7    277      1
## 3 01/07/2003     99     29    335     26 18/05/2005     57     51    342     42
## 4 01/10/2003    188      8    330     12 18/08/2005    145     51     41     25
## 5 01/01/2002    281     10    347     26 19/11/2003    236     58    345     29
## 6 01/04/2002     11     50     51     47 17/02/2004    328     21     39      5
## 7 01/07/2002     99     44    112     20 18/05/2004     58      3     97      9
## 8 01/10/2002    188     23    171      5 18/08/2004    146      3    155     13
ls1 = pi *(mydata2$ls1deg + mydata2$ls1min/60)/180.
ls2 = pi *(mydata2$ls2deg + mydata2$ls2min/60)/180.
lm1 = pi *(mydata2$lm1deg + mydata2$lm1min/60)/180.
lm2 = pi *(mydata2$lm2deg + mydata2$lm2min/60)/180.

Build the list of positions of Mars on its orbit:

xm
## [1] -0.80266
ym
## [1] -0.8599597
getmars <-function(i){
  
  am=-1.0100
  bm=-0.003108
  cm=-0.2063
  dm=-0.006017
  fm= 1.2065
  
dmi=10000
theta1 = ls1[i]-lm0+pi
aa=a*(tan(theta1))^2 + b*tan(theta1)-1
bb=c+d*tan(theta1)
cc=f
delta=bb*bb-4*aa*cc
x1a=(-bb+sqrt(delta))/(2*aa)
y1a=x1a*tan(theta1)
x1b=(-bb-sqrt(delta))/(2*aa)
y1b=x1b*tan(theta1)
theta2 = ls2[i]-lm0+pi
aa=a*(tan(theta2))^2 + b*tan(theta2)-1
bb=c+d*tan(theta2)
cc=f
delta=bb*bb-4*aa*cc
x2a=(-bb+sqrt(delta))/(2*aa)
y2a=x2a*tan(theta2)
x2b=(-bb-sqrt(delta))/(2*aa)
y2b=x2b*tan(theta2)




#ax x + ay y = az
#bx x + by y = bz 

# case a a
ax=tan(lm1[i]-lm0)
ay=-1
az=-y1a+x1a*tan(lm1[i]-lm0)
bx=tan(lm2[i]-lm0)
by=-1
bz=-y2a+x2a*tan(lm2[i]-lm0)
m=matrix(c(ax,ay,bx,by),nrow=2,byrow=TRUE)
mi=solve(m)
rhs=matrix(c(az,bz),nrow=2,byrow=TRUE)
v=mi%*%rhs

xx=v[1]
yy=v[2]
dse=-xx^2+am*yy^2+bm*xx*yy+cm*xx+dm*yy+fm
  vsol=v
  dmi=dse

# case a b
ax=tan(lm1[i]-lm0)
ay=-1
az=-y1a+x1a*tan(lm1[i]-lm0)
bx=tan(lm2[i]-lm0)
by=-1
bz=-y2b+x2b*tan(lm2[i]-lm0)
m=matrix(c(ax,ay,bx,by),nrow=2,byrow=TRUE)
mi=solve(m)
rhs=matrix(c(az,bz),nrow=2,byrow=TRUE)
v=mi%*%rhs

xx=v[1]
yy=v[2]
dse=-xx^2+am*yy^2+bm*xx*yy+cm*xx+dm*yy+fm
if(dse <= dmi){
  vsol=v
  dmi=dse
}

# case b a
ax=tan(lm1[i]-lm0)
ay=-1
az=-y1b+x1b*tan(lm1[i]-lm0)
bx=tan(lm2[i]-lm0)
by=-1
bz=-y2a+x2a*tan(lm2[i]-lm0)
m=matrix(c(ax,ay,bx,by),nrow=2,byrow=TRUE)
mi=solve(m)
rhs=matrix(c(az,bz),nrow=2,byrow=TRUE)
v=mi%*%rhs

xx=v[1]
yy=v[2]
dse=-xx^2+am*yy^2+bm*xx*yy+cm*xx+dm*yy+fm
if(dse <= dmi){
  vsol=v
  dmi=dse
}


# case b b
ax=tan(lm1[i]-lm0)
ay=-1
az=-y1b+x1b*tan(lm1[i]-lm0)
bx=tan(lm2[i]-lm0)
by=-1
bz=-y2b+x2b*tan(lm2[i]-lm0)
m=matrix(c(ax,ay,bx,by),nrow=2,byrow=TRUE)
mi=solve(m)
rhs=matrix(c(az,bz),nrow=2,byrow=TRUE)
v=mi%*%rhs

xx=v[1]
yy=v[2]
dse=-xx^2+am*yy^2+bm*xx*yy+cm*xx+dm*yy+fm
if(dse <= dmi){
  vsol=v
  dmi=dse
}

print(vsol)

return(vsol)
}
gett1(1,1)
## [1] 1.000000 1.000000 4.902630 4.014839
## [1] -0.4176677  0.5762791
gett1(1,-1)
## [1]  1.000000 -1.000000  4.902630  4.014839
## [1]  0.4323007 -0.5964691
gett1(4,1)
## [1] 4.000000 1.000000 3.283546 5.763077
## [1] -0.6047123 -0.3952103
gett1(4,-1)
## [1]  4.000000 -1.000000  3.283546  5.763077
## [1] 0.6067002 0.3965096
plot(xyDirect[,1],xyDirect[,2],type='l',
col='cyan',xlim=c(-1,1),ylim=c(-1,1),asp=1);par(new=TRUE)
lines(c(0,1),c(0,0))


vv1=gett1(1,1)
## [1] 1.000000 1.000000 4.902630 4.014839
lines(c(0,vv1[1]),c(0,vv1[2]),col='red')

t=lm1[1]-lm0
uu=c(cos(t),sin(t))
m1=vv1
m2=m1+2*uu

lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')

vv1=gett1(1,-1)
## [1]  1.000000 -1.000000  4.902630  4.014839
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')

vv2=gett2(1,1)
## [1] 1.000000 1.000000 4.130613 3.751585
lines(c(0,vv2[1]),c(0,vv2[2]),col='pink')

vv2=gett2(1,-1)
## [1]  1.000000 -1.000000  4.130613  3.751585
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

t=lm2[1]-lm0
uu=c(cos(t),sin(t))

m1=vv2
m2=m1+2*uu

lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')

i=2
plot(xyDirect[,1],xyDirect[,2],type='l',
col='cyan',xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),asp=1);par(new=TRUE)
lines(c(0,1),c(0,0))


vv1=gett1(i,1)
## [1] 2.0000000 1.0000000 0.2024582 5.0189852
lines(c(0,vv1[1]),c(0,vv1[2]),col='red')


vv1=gett1(i,-1)
## [1]  2.0000000 -1.0000000  0.2024582  5.0189852
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')


t=lm1[1]-lm0
uu=c(cos(t),sin(t))
m1=vv1
m2=m1+2*uu

lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')



vv2=gett2(i,1)
## [1] 2.000000 1.000000 5.726716 4.834853
lines(c(0,vv2[1]),c(0,vv2[2]),col='pink')

vv2=gett2(i,-1)
## [1]  2.000000 -1.000000  5.726716  4.834853
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

t=lm2[1]-lm0
uu=c(cos(t),sin(t))

m1=vv2
m2=m1+2*uu

lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')

plot(xyDirect[,1],xyDirect[,2],type='l',
col='cyan',xlim=c(-2,2),ylim=c(-2,2));par(new=TRUE)
lines(c(0,1),c(0,0))
vv=gett1(1,1)
## [1] 1.000000 1.000000 4.902630 4.014839
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(1,-1)
## [1]  1.000000 -1.000000  4.902630  4.014839
lines(c(0,vv[1]),c(0,vv[2]),col='blue')

vv=gett1(2,1)
## [1] 2.0000000 1.0000000 0.2024582 5.0189852
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(2,-1)
## [1]  2.0000000 -1.0000000  0.2024582  5.0189852
lines(c(0,vv[1]),c(0,vv[2]),col='blue')


vv=gett1(3,1)
## [1] 3.000000 1.000000 1.736312 5.854416
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(3,-1)
## [1]  3.000000 -1.000000  1.736312  5.854416
lines(c(0,vv[1]),c(0,vv[2]),col='blue')


vv=gett1(4,1)
## [1] 4.000000 1.000000 3.283546 5.763077
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(4,-1)
## [1]  4.000000 -1.000000  3.283546  5.763077
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(5,1)
## [1] 5.000000 1.000000 4.907284 6.063856
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(5,-1)
## [1]  5.000000 -1.000000  4.907284  6.063856
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(6,1)
## [1] 6.0000000 1.0000000 0.2065306 0.9037897
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(6,-1)
## [1]  6.0000000 -1.0000000  0.2065306  0.9037897
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(7,1)
## [1] 7.000000 1.000000 1.740675 1.960587
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(7,-1)
## [1]  7.000000 -1.000000  1.740675  1.960587
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(8,1)
## [1] 8.000000 1.000000 3.287909 2.985967
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(8,-1)
## [1]  8.000000 -1.000000  3.287909  2.985967
lines(c(0,vv[1]),c(0,vv[2]),col='blue')


vv=gett2(1,1)
## [1] 1.000000 1.000000 4.130613 3.751585
lines(c(0,vv[1]),c(0,vv[2]),col='pink')
vv=gett2(1,-1)
## [1]  1.000000 -1.000000  4.130613  3.751585
lines(c(0,vv[1]),c(0,vv[2]),col='green')

#lines(gett1(2,1),gett1(2,-1))
#lines(gett1(4,1),gett1(4,-1))
getmars <-function(i){
  
  am=-1.0100
  bm=-0.003108
  cm=-0.2063
  dm=-0.006017
  fm= 1.2065
  
  a=-0.9996
  b=-0.0002594
  c= 0.01495
  d=-0.01926
  f= 0.5242
  
dmi=10000
theta1 = ls1[i]-lm0+pi
aa=a*(tan(theta1))^2 + b*tan(theta1)-1
bb=c+d*tan(theta1)
cc=f
delta=bb*bb-4*aa*cc
x1a=(-bb+sqrt(delta))/(2*aa)
y1a=x1a*tan(theta1)
dse=sqrt(x1a^2+y1a^2)

dto=abs(0.72-dse)

if(dto <= dmi){
  x1=x1a
  y1=y1a
  dmi=dto
}

x1b=(-bb-sqrt(delta))/(2*aa)
y1b=x1b*tan(theta1)

dse=sqrt(x1b^2+y1b^2)

dto=abs(0.72-dse)

if(dto <= dmi){
  x1=x1b
  y1=y1b
  dmi=dto
}


dmi=10000
theta2 = ls2[i]-lm0+pi
aa=a*(tan(theta2))^2 + b*tan(theta2)-1
bb=c+d*tan(theta2)
cc=f
delta=bb*bb-4*aa*cc
x2a=(-bb+sqrt(delta))/(2*aa)
y2a=x2a*tan(theta2)

dse=sqrt(x2a^2+y2a^2)

dto=abs(0.72-dse)

if(dto <= dmi){
  x2=x2a
  y2=y2a
  dmi=dto
}

x2b=(-bb-sqrt(delta))/(2*aa)
y2b=x2b*tan(theta2)
dse=sqrt(x2b^2+y2b^2)

dto=abs(0.72-dse)

if(dto <= dmi){
  x2=x2b
  y2=y2b
  dmi=dto
}

#ax x + ay y = az
#bx x + by y = bz 

ax=tan(lm1[i]-lm0)
ay=-1
az=-y1+x1*tan(lm1[i]-lm0)
bx=tan(lm2[i]-lm0)
by=-1
bz=-y2+x2*tan(lm2[i]-lm0)
m=matrix(c(ax,ay,bx,by),nrow=2,byrow=TRUE)
mi=solve(m)
rhs=matrix(c(az,bz),nrow=2,byrow=TRUE)
vsol=mi%*%rhs

return(vsol)
}
xm=c()
ym=c()

for(i in 1:8){

vv=getmars(i)
xm=c(xm,vv[1])
ym=c(ym,vv[2])
  
}

xy <- matrix(c(xm,ym),8,2,byrow=FALSE)
xy
##             [,1]       [,2]
## [1,] -0.79954725 -0.8545204
## [2,]  0.02776072 -1.0936782
## [3,] -0.84602172  0.5831269
## [4,] -0.91662227 -0.3695864
## [5,]  0.63911007  0.8088850
## [6,] -0.22216479  1.0894548
## [7,]  0.95172769 -0.6650414
## [8,]  1.18685046  0.1191178
plot(xy)

xym <- matrix(c(xm,ym),8,2,byrow=FALSE)
xym
##             [,1]       [,2]
## [1,] -0.79954725 -0.8545204
## [2,]  0.02776072 -1.0936782
## [3,] -0.84602172  0.5831269
## [4,] -0.91662227 -0.3695864
## [5,]  0.63911007  0.8088850
## [6,] -0.22216479  1.0894548
## [7,]  0.95172769 -0.6650414
## [8,]  1.18685046  0.1191178
plot(xym)

sqrt(xm^2+ym^2)
## [1] 1.1702483 1.0940304 1.0275163 0.9883271 1.0309009 1.1118763 1.1610623
## [8] 1.1928131
ellipDirectm <- EllipseDirectFit(xy)
ellipDirectGm <- AtoG(ellipDirect)$ParG
xyDirectm<-calculateEllipse(ellipDirectGm[1], ellipDirectGm[2], ellipDirectGm[3],
ellipDirectGm[4], 180/pi*ellipDirectGm[5])
plot(xyDirectm[,1],xyDirectm[,2],type='l',
col='cyan');par(new=TRUE)
lines(xym)

ellipDirectm
##             [,1]
## [1,] -0.53700753
## [2,] -0.00193463
## [3,] -0.53695671
## [4,]  0.05595009
## [5,] -0.03600174
## [6,]  0.64720176
print('a')
## [1] "a"
bm=-ellipDirectm[2]/ellipDirectm[1]
am=-ellipDirectm[3]/ellipDirectm[1]
cm=-ellipDirectm[4]/ellipDirectm[1]
dm=-ellipDirectm[5]/ellipDirectm[1]
fm=-ellipDirectm[6]/ellipDirectm[1]

am
## [1] -0.9999054
bm
## [1] -0.003602612
cm
## [1] 0.1041886
dm
## [1] -0.06704141
fm
## [1] 1.205201
-xm^2+am*ym^2+bm*xm*ym+cm*xm+dm*ym+fm
## [1] -0.192688256  0.084734619  0.023980724  0.156478469  0.153002624
## [6] -0.126269691  0.003201501 -0.102439990
XY <- matrix(c(xm,ym),8,2,byrow=FALSE)
XY
##             [,1]       [,2]
## [1,] -0.79954725 -0.8545204
## [2,]  0.02776072 -1.0936782
## [3,] -0.84602172  0.5831269
## [4,] -0.91662227 -0.3695864
## [5,]  0.63911007  0.8088850
## [6,] -0.22216479  1.0894548
## [7,]  0.95172769 -0.6650414
## [8,]  1.18685046  0.1191178
plot(XY)

xy=XY
XY <- matrix(c(xm,ym),8,2,byrow=FALSE)
XY
##             [,1]       [,2]
## [1,] -0.79954725 -0.8545204
## [2,]  0.02776072 -1.0936782
## [3,] -0.84602172  0.5831269
## [4,] -0.91662227 -0.3695864
## [5,]  0.63911007  0.8088850
## [6,] -0.22216479  1.0894548
## [7,]  0.95172769 -0.6650414
## [8,]  1.18685046  0.1191178
plot(XY)

xy=XY


#EllipseDirectFit(XY)

ellipDirect <- EllipseDirectFit(xy)
ellipDirectG <- AtoG(ellipDirect)$ParG
xyDirect<-calculateEllipse(ellipDirectG[1], ellipDirectG[2], ellipDirectG[3],
ellipDirectG[4], 180/pi*ellipDirectG[5])
plot(xyDirect[,1],xyDirect[,2],type='l',
col='cyan');par(new=TRUE)
lines(XY)

ellipDirect
##             [,1]
## [1,] -0.53700753
## [2,] -0.00193463
## [3,] -0.53695671
## [4,]  0.05595009
## [5,] -0.03600174
## [6,]  0.64720176
print('a')
## [1] "a"
b=-ellipDirect[2]/ellipDirect[1]
a=-ellipDirect[3]/ellipDirect[1]
c=-ellipDirect[4]/ellipDirect[1]
d=-ellipDirect[5]/ellipDirect[1]
f=-ellipDirect[6]/ellipDirect[1]

a
## [1] -0.9999054
b
## [1] -0.003602612
c
## [1] 0.1041886
d
## [1] -0.06704141
f
## [1] 1.205201

The original work of Kepler

https://www.keplersdiscovery.com/Earth.html

Principle of the method

Sidereal period of Mars

Kepler believed in the system of Copernicus. But how long did it take for Mars to complete a revolution around the Sun? If we use far away stars to build a rigid frame, how much time will it take for Mars to come back to the same position in this frame?

The problem is: we do not sit on a fixed star, but on a moving planet (Earth) and observe Mars from this moving platform. Hence we must somehow find a way to remove this motion to get access to the sidereal period of Mars, i.e. the number of terrestrial days after which Mars will have completed a revolution around the Sun, seen from stars. This period is obtained by making the following observation.

On one hand, we know how much time in takes (one year or 365.25 days) for Earth to complete its revolution around the Sun (by observing the apparent motion of the Sun in the sky and its recurrence). This is an information involving the relative motion of the Earth with respect to the Sun. But we also can observe from Earth the time it takes for an opposition of Mars to take place. Mars in in opposition when the Earth, the Sun and Mars are aligned, with the Earth in between. At oppositions, the planet rises when the sun sets.

The interval of time between oppositions is called the synodic period, and its value for Mars is \(S_M=780d\).

Now, express this period in terrestrial and martian years.

It takes \((1+\alpha)\) martian years to complete a synodic revolution:

\[S_M=(1+\alpha)T_M\] This is equivalent to saying that the planet Mars will have covered one revolution and a fraction \(\alpha\) of revolution around the Sun. If we start with an alignment STM, to arrive at another alignment, the Earth will have had to cover two revolutions and the same fraction, or

\[S_M=(2+\alpha)T_T\] This translates, by elimination of \(\alpha\), into:

\[\frac{S_M}{T_M} -1 = \frac{S_M}{T_T} - 2 \]

Or

\[\frac{1}{T_M} = \frac{1}{T_T} - \frac{1}{S_M}\]

Note that the signs would be reversed for inferior planets (Venus and Mercury), but here the sign is correct: Mars is a superior planet.

Because \(T_T=365.25d\) and \(S_M=780d\), we arrive at the following value for the sidereal period of Mars:

\[T_M=687d\]

Triangulation of positions

The genial idea of Kepler was to use the sidereal period of Mars to triangulate Mars’s and Earth’s orbits by ‘stroboscopy’.

We now know that Mars returns at the same position every 687 days. This means that we can view the same position (of Mars) from Earth from different perspectives. Which means that we can reconstruct where Mars is (same but unknown position) from different line of sights as long as these correspond to observations taken at intervals of \(T_M=687\) days.

Geometric construction by hand

I invite you to try the geometric construction by hand, as explained in this

PHYS1401 Mars Lab

And this for two reasons: -it gives a good intuition about the method devised by Kepler -it is not very accurate and shows the limits of hand-drawn constructions

The errors involved in drawing each position of Mars are likely to lead to a sub-optimal use of the data. It would be interesting to quantify them. But what is the point of loosing the benefit of laborious and precise observations (longitudes were recorded by Tycho Brahe to a precision of within one minute of arc) if we should afterwards use a slopy method?

The heliocentric longitude of Earth is the position of the Earth relative to the sun, and the geocentric longitude of Mars is the position of Mars relative to the Earth. Use the following table of observations to construct the orbit of the planet Mars:

Note how these ‘observations’ are only within 1 degree!

DATE HELIOCENTRIC LONGITUDE OF EARTH GEOCENTRIC LONGITUDE OF MARS 1 March 21, 1931 180° 118° February 5, 1933 136° 168° 2 April 20, 1933 210° 151° March 8, 1935 167° 204° 3 May 26, 1935 245° 187° April 12, 1937 202° 246° 4 September 26, 1939 3° 301° August 13, 1941 320° 20° 5 November 22, 1941 60° 12° October 10, 1943 17° 80° 6 January 21, 1944 121° 66° December 8, 1945 76° 123° 7 March 19, 1946 178° 108° February 4, 1948 135° 153° 8 April 4, 1948 195° 138° February 20, 1950 152° 191°

Determination of geocentric ecliptic longitudes

By calculation

Fitzpatrick (New Almagest)

The geometry of the Keplerian orbit is shown in the following figure:

The orbital eccentricity is \(e\). The angle \(SCQ=E\) is called the ecliptic anomaly. The angle \(RSP=T\) is the true anomaly.

Kepler defined the mean anomaly, \(M\), which is an angle that increases at constant speed and is in one-to-one relation with the position of the satellite on its orbit.

It is defined as:

\[M=\frac{2\pi}{\tau}(t-t^{\ast})\]

The period of this motion is the period of the satellite, \(\tau\), i.e. the time that the satellites takes to complete one revolution.

\[M=E-e\sin E\]

\[\begin{eqnarray} \frac{r}{a} &=& 1 -e\,\cos M + e^2\,\sin^2 M,\\[0.5ex] T &=& M + 2\,e\,\sin M + (5/4)\,e^2\,\sin 2M. \end{eqnarray}\]

difference between https://stackoverflow.com/questions/31843662/what-is-the-difference-between-cat-and-print https://stat.ethz.ch/pipermail/r-help/2004-September/056974.html

paste('u = ',4)
## [1] "u =  4"
mydata=read.csv("https://mathstat.dal.ca/~fullsack/DATA/pldata1.csv")
summary(mydata)
##    PeriodDays       AbsMag1          AbsMag2      
##  Min.   : 2.00   Min.   :-6.400   Min.   :-4.800  
##  1st Qu.: 6.25   1st Qu.:-5.725   1st Qu.:-4.050  
##  Median :12.50   Median :-4.300   Median :-2.700  
##  Mean   :28.28   Mean   :-4.500   Mean   :-2.894  
##  3rd Qu.:47.50   3rd Qu.:-3.550   3rd Qu.:-1.950  
##  Max.   :90.00   Max.   :-2.300   Max.   :-0.800
head(mydata)
##   PeriodDays AbsMag1 AbsMag2
## 1          2    -2.3    -0.8
## 2          3    -2.8    -1.2
## 3          4    -3.1    -1.5
## 4          5    -3.3    -1.7
## 5          6    -3.5    -1.9
## 6          7    -3.7    -2.1

By measurements