---
title: "Kepler and Mars"
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---


# Outline

In this activity, we shall apply the principles of regression to the problem of accurate reconstruction of the orbits of planets. We will in fact almost repeat exactly the method invented by the german astronomer Jonathan Kepler to reconstruct the orbit of Earth and of Mars around the Sun, based on astronomical observations of Tycho Brahe (the 'data').

This work took Kepler eight years to accomplish. Hopefully, this activity should take you less time.

The plan of this activity is as follows. We will first use Kepler's approach to reconstruct the orbit of the Earth. We will then attempt to reconstruct the orbit of Mars but will meet some computational problem.
In the next section, we will elegantly overcome this problem by recasting it into a more proper form  - using adequate intersection primitives/operations. This shall be entirely satisfactory and shed a bright light on Kepler's method.


Please note that I originally devised this activity for a worskop in astronomy destined to the students of the Alliance Francaise Halifax school.

(To those who might be tempted: please do not vampirize or sell this activity online...) 


## Learning objectives

This activity has so many learning objectives that I renounce to make an exhaustive list. Suffice it to mention that we are going to use state-of-the-art R packages and Rstudio to explore a most exemplary episod of scientific revolution. Using a least-square methods to fit conics in this context should relieve you from the deeply boring and more conventional exercises on linear regression (see my lectures for those who like this stuff). My philosophy about statistics is that it is often the least interesting part of science. It should be clear that statistical analysis is only useful when a sound framework exists for collecting observations. The cleverness of Kepler is not related to stats (he did not use any) but to his ingenious selection and a priori understanding of the data of Tycho. However regression and least square methodshave the merit of providing a universal method for reducing errors and can tackle dirtier and more subtle -if lss fundamental- problems (for exemple the tracking and forecasting of satellite or comet orbit in the presence of friction forces or ice abrasion).


# First attempt



# Second attempt

```{r}
#install.packages('PlaneGeometry')
library(PlaneGeometry)
C   =c(0.1,-0.2)
S   =c(0,0)
gam =c(0,-1)


#Line$new(A, B, extendA = TRUE, extendB = TRUE)
ell <- Ellipse$new(c(1,1), 5, 1, 30)
line <- Line$new(c(2,-2), c(0,4))
( Is <- intersectionEllipseLine(ell, line) )
ell$includes(Is$I1); ell$includes(Is$I2)

circ <- Circle$new(c(1,1), 2)
line <- Line$new(c(2,-2), c(1,2), FALSE, FALSE)
intersectionCircleLine(circ, line)
intersectionCircleLine(circ, line, strict = TRUE)

#intersectionLineLine(line1, line2, strict = FALSE)
```

### define the terrestrial orbit

let us learn how to do plots...

```{r}

#plotrix
#draw.circle(x,y,radius,nv=100,border=NULL,col=NA,lty=1,density=NULL,
#  angle=45,lwd=1)



# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(0,2.5), ylim = c(0,2.5),
xlab = NA, ylab = NA)
grid()
# draw a triangle
t <- Triangle$new(c(0,0), c(1,0), c(0.5,sqrt(3)/2))
draw(t, col = "blue", lwd = 2)
draw(t$rotate(90, t$C), col = "green", lwd = 2)
# draw a circle
circ <- t$incircle()
draw(circ, col = "orange", border = "brown", lwd = 2)
# draw an ellipse
S <- Scaling$new(circ$center, direction = c(2,1), scale = 2)
draw(S$scaleCircle(circ), border = "grey", lwd = 2)
# draw a line
l <- Line$new(c(1,1), c(1.5,1.5), FALSE, TRUE)
draw(l, col = "red", lwd = 2)
perp <- l$perpendicular(c(2,1))
draw(perp, col = "yellow", lwd = 2)


```


We are going to look at the earth orbit from above
kepler new that the motion of earth around the sun was in a plane- effectively the apparent motion of the sun seen from Earth is in a plane, the plane of the ecliptic, so called because all eclipses of the sun arise in or near this plane.
This is because all planets orbit around the sun in planes close to the ecliptic plane, the plane that contains the earth orbit.



```{r}
e=0.01653
a=1/e

# I need to explain this angle
# this is the ecliptic longitude of the perihelion
@ value will vary with epochs, not the same today and at the time of Kepler
omega=100+19/60
omega=pi*omega/180

om=omega+pi
C   =c(0.1,-0.2)
S   =c(0,0)
gamnor =c(-gam[2],gam[1])

ep=0.3
dir=S+cos(om)*gam+sin(om)*gamnor

# center of circular earth orbit, close to the Sun and on the line of apsides
C=S+ep*(cos(om)*gam+sin(om)*gamnor)
# aphelion 
A=S+(1+ep)*(cos(om)*gam+sin(om)*gamnor)
# perihelion
P=S-(1-ep)*(cos(om)*gam+sin(om)*gamnor)

#plotrix
#draw.circle(x,y,radius,nv=100,border=NULL,col=NA,lty=1,density=NULL,
#  angle=45,lwd=1)

# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(-1,1), ylim = c(-1,1),
xlab = NA, ylab = NA,main="Earth orbit with exagerated eccentricity")
#grid()
# draw a circle
#circ <- t$incircle()
#draw(circ, col = "orange", border = "brown", lwd = 2)

# draw a line
l <- Line$new(S, gam, FALSE, TRUE)
draw(l, col = "red", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

l <- Line$new(C, A, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(C, P, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

circ <- Circle$new(C, 1)
draw(circ, border = "brown", lwd = 2)

```



```{r}
e=0.01653
a=1/e
omega=100+19/60
omega=pi*omega/180

om=omega+pi
C   =c(0.1,-0.2)
S   =c(0,0)
gamnor =c(-gam[2],gam[1])

ep=e
dir=S+cos(om)*gam+sin(om)*gamnor
C=S+ep*(cos(om)*gam+sin(om)*gamnor)

A=S+(1+ep)*(cos(om)*gam+sin(om)*gamnor)
P=S-(1-ep)*(cos(om)*gam+sin(om)*gamnor)


  

#plotrix
#draw.circle(x,y,radius,nv=100,border=NULL,col=NA,lty=1,density=NULL,
#  angle=45,lwd=1)

# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(-1,1), ylim = c(-1,1),
xlab = NA, ylab = NA,main="Earth orbit with true eccentricity")
#grid()
# draw a circle
#circ <- t$incircle()
#draw(circ, col = "orange", border = "brown", lwd = 2)

# draw a line
l <- Line$new(S, gam, FALSE, TRUE)
draw(l, col = "red", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

l <- Line$new(C, A, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(C, P, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

circT <- Circle$new(C, 1)
draw(circT, border = "brown", lwd = 2)



```


This may look like a small perturbation compared to a cercle centered on the sun, but it is a critical one. No accurate prediction of astronomic events can be made if you neglect this little blue segment - i.e. if you zero the Earth eccentricity.
Note however that this model of Kepler does not quite represent the truth: the Earth orbit is, like the orbit of Mars, an ellipse, not a circle. But the trace of the true elliptic orbit and of this off-centered circular orbit are very close. The error made is of the order of th radius of the earth.

orbit of Mars and Earth

<img src="https://www.eg.bucknell.edu/physics/astronomy/astr101/calpages/wom4.gif">





Now let us prove that the intersection primitives work



```{r}

lont1a=159+38/60
lont1a=pi*lont1a/180

dirt1a=S+cos(lont1a)*gam+sin(lont1a)*gamnor

#plotrix
#draw.circle(x,y,radius,nv=100,border=NULL,col=NA,lty=1,density=NULL,
#  angle=45,lwd=1)

# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(-1,1), ylim = c(-1,1),
xlab = NA, ylab = NA,main="Earth orbit with true eccentricity")
#grid()
# draw a circle
#circ <- t$incircle()
#draw(circ, col = "orange", border = "brown", lwd = 2)

# draw a line
l <- Line$new(S, gam, FALSE, TRUE)
draw(l, col = "red", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

l <- Line$new(C, A, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(C, P, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

circT <- Circle$new(C, 1)
draw(circT, border = "brown", lwd = 2)
# draw a line
l <- Line$new(S, dirt1a, FALSE, TRUE)
draw(l, col = "orange", lwd = 2)


t1a=intersectionCircleLine(circT, l, strict = TRUE)

lt1a <- Line$new(P, t1a, FALSE, FALSE)
draw(lt1a, col = "orange", lwd = 2)

```

Now let us build the first stage of our scaffold. Ww will reconstruct the two positions of Earth on iits orbit for the first pair of data


```{r}

lont1a=159+38/60
lont1a=pi*lont1a/180

dirt1a=S+cos(lont1a)*gam+sin(lont1a)*gamnor


lont1b=115+35/60
lont1b=pi*lont1b/180

dirt1b=S+cos(lont1b)*gam+sin(lont1b)*gamnor


# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(-1,1), ylim = c(-1,1),
xlab = NA, ylab = NA,main="Scaffolding: step 1")
#grid()
# draw a circle
#circ <- t$incircle()
#draw(circ, col = "orange", border = "brown", lwd = 2)

# draw a line
l <- Line$new(S, gam, FALSE, TRUE)
draw(l, col = "red", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

l <- Line$new(C, A, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(C, P, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

circT <- Circle$new(C, 1)
draw(circT, border = "brown", lwd = 2)
# draw a line
l1a <- Line$new(S, dirt1a, FALSE, TRUE)
#draw(l, col = "orange", lwd = 2)


t1a=intersectionCircleLine(circT, l1a, strict = TRUE)

lt1a <- Line$new(S, t1a, FALSE, FALSE)
draw(lt1a, col = "orange", lwd = 2)

l1b <- Line$new(S, dirt1b, FALSE, TRUE)
t1b=intersectionCircleLine(circT, l1b, strict = TRUE)

lt1b <- Line$new(S, t1b, FALSE, FALSE)
draw(lt1b, col = "orange", lwd = 2)

```

Now let us build the next floor of our scaffold:


```{r}

lonm1a=135+20/60
lonm1a=pi*lonm1a/180

dirm1a=t1a+cos(lonm1a)*gam+sin(lonm1a)*gamnor


lonm1b=182+13/60
lonm1b=pi*lonm1b/180

dirm1b=t1b+cos(lonm1b)*gam+sin(lonm1b)*gamnor


# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(-2,2), ylim = c(-2,2),
xlab = NA, ylab = NA,main="Scaffolding: step 2")
#grid()
# draw a circle
#circ <- t$incircle()
#draw(circ, col = "orange", border = "brown", lwd = 2)

# draw a line
l <- Line$new(S, gam, FALSE, TRUE)
draw(l, col = "red", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

l <- Line$new(C, A, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(C, P, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

circT <- Circle$new(C, 1)
draw(circT, border = "brown", lwd = 2)
# draw a line
l1a <- Line$new(S, dirt1a, FALSE, TRUE)
#draw(l, col = "orange", lwd = 2)


t1a=intersectionCircleLine(circT, l1a, strict = TRUE)

lt1a <- Line$new(S, t1a, FALSE, FALSE)
draw(lt1a, col = "orange", lwd = 2)

l1b <- Line$new(S, dirt1b, FALSE, TRUE)
t1b=intersectionCircleLine(circT, l1b, strict = TRUE)

lt1b <- Line$new(S, t1b, FALSE, FALSE)
draw(lt1b, col = "orange", lwd = 2)


lm1a <- Line$new(t1a, dirm1a, FALSE, TRUE)
draw(lm1a, col = "blue", lwd = 2)



```





```{r}

plotbas <-function(){
# open new plot window
plot(0, 0, type="n", asp = 1, xlim = c(-2,2), ylim = c(-2,2),
xlab = NA, ylab = NA,main="Scaffolding: step 2")
#grid()
# draw a circle
#circ <- t$incircle()
#draw(circ, col = "orange", border = "brown", lwd = 2)

# draw a line
l <- Line$new(S, gam, FALSE, TRUE)
draw(l, col = "red", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

l <- Line$new(C, A, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(C, P, FALSE, FALSE)
draw(l, col = "green", lwd = 2)
l <- Line$new(S, C, FALSE, FALSE)
draw(l, col = "blue", lwd = 2)

circT <- Circle$new(C, 1)
draw(circT, border = "brown", lwd = 2)

}

plotbas()

lm1a <- Line$new(t1a, dirm1a, FALSE, TRUE)
draw(lm1a, col = "blue", lwd = 2)



```





```{r}

plotbas()

lm1a <- Line$new(t1a, dirm1a, FALSE, TRUE)
draw(lm1a, col = "blue", lwd = 2)

lm1b <- Line$new(t1b, dirm1b, FALSE, TRUE)
draw(lm1b, col = "blue", lwd = 2)

```


Fantastic we just pinned Mars in space !

Let us use the intersection primitives to get the reconstructed position of Mars on its orbit:


```{r}
plotbas()


m1=intersectionLineLine(lm1a, lm1b, strict = TRUE)

lim1a <- Line$new(t1a, m1, FALSE, FALSE)
draw(lim1a, col = "blue", lwd = 2)

lim1b <- Line$new(t1b, m1, FALSE, FALSE)
draw(lim1b, col = "blue", lwd = 2)

```



## Martian Ode to Joy


To express our joy on such a magnificent discovery, we can imagine Mars, god of war, or maybe how the Rudolph II, Emperor of the Bohemian Empire and protector of Kepler at the court of Prague, must have felt upon hearing so magnificent a discovery.


<img src="">

<img width=500 src="https://upload.wikimedia.org/wikipedia/commons/thumb/b/b5/Vertumnus_%C3%A5rstidernas_gud_m%C3%A5lad_av_Giuseppe_Arcimboldo_1591_-_Skoklosters_slott_-_91503.tiff/lossy-page1-826px-Vertumnus_%C3%A5rstidernas_gud_m%C3%A5lad_av_Giuseppe_Arcimboldo_1591_-_Skoklosters_slott_-_91503.tiff.jpg">


portrait mythologique de henri 4 en dieu mars
https://www.jstor.org/stable/43445017?seq=1#metadata_info_tab_contents

link

<img width=500 src="https://upload.wikimedia.org/wikipedia/commons/a/af/Mus%C3%A9e_national_du_Ch%C3%A2teau_de_Pau_-_Portait_d%27Henri_IV_en_Mars_-_Ambroise_Dubois_P_81_20_1.jpg">

<a href="43445017.pdf">henri 4 en mars </a>

Représentation allégorique d'Henri IV rex imperator
MARIE-FRANCE WAGNER
C'est à partir de deux images que nous étudions la représentation d'Henri IV rex imperator. La première, gravée en 1602 et réutilisée quelque sept années plus tard, est doublement allégorique; elle met en scène Henri IV à la fois Hercule et Alexandre. La description du premier arc de triomphe de l'Entrée en Avignon de 1600 peut lui servir de support scriptural. La seconde est un portrait mythologique d'Henri IV en Mars, qui foule à ses pieds ses armes devenues inutiles. Ces images symboliques s'efforcent de présenter le roi converti, puis sacré, dans sa fonction de pacificateur (vainqueur des guerres civiles politico-religieuses) et de restaurateur de la monarchie. Elles représentent la grandeur 


la representation en chef de guerre de rodophe ii ds la peint suiv n'est pas ss rappeller le dieu mars et
la fig all de henri 4

Portrait de l'empereur Rodolphe II
1604
Aegidius Sadeler II
Hans van Aachen Hans von Aachen (Cologne, 1552 – Prague, 4 mars 1615 ) est un peintre maniériste allemand

<img width=500 src="https://collections.geneve.ch/mah/sites/default/files/styles/scale_w1148/public/image/2019-11/265941.jpeg?itok=tx6WiNpD">



the planets venus and mars impregnate  the collective imagination

The often strongly eroticized portraits of Mars and Venus form the basis of a fundamental equilibrium of forces, a form of equivalent of the oriental Ying and Yang for the occidental civilization of the classical age. The symbolism of mannerism

Mannerism is characterized by
The mannerist representations of Venus and Mars are far remote from the idealized aetheral esthetic of the renaissance paintings of Boticelli.
Venus and Mars (or Mars and Venus) is a panel painting of about 1485 by the Italian Renaissance painter Sandro Botticelli


leaving indelible marks on the collective imagination.



Eros et poesia, la peinture a la cour
Mars, Venus and Amor

Bartholomäus Spranger

<img src="">

https://musee-magnin.fr/sites/musee-magnin.fr/files/dp_praga_magica.pdf



```{r}

plotbas()

lm1a <- Line$new(t1a, dirm1a, FALSE, TRUE)
draw(lm1a, col = "blue", lwd = 2)

lm1b <- Line$new(t1b, dirm1b, FALSE, TRUE)
draw(lm1b, col = "blue", lwd = 2)

```




# A primer on the sky seen from Earth

### The ecliptic plane


<img src="">


### The vernal point

https://sites.google.com/site/astrobazar/le-pied-a-l-etrier/6-sky-high-way/image004.jpg?attredirects=0

<img src="">


### Retrograde motion of planets

In the cosmological model of Aristotle, all celestial bodies rotate around the Earth, which is the center of the universe. However, some observational facts do not support this notion. For example, planets seem to occasionally travel backwards in the sky.

Here is a superposition of photographs of Mars seen from Earth showing a retogradation loop.

<img src="http://www.erasmatazz.com/_Media/retromars_med_hr.jpeg">


The roman/egyptian astronomer Claudius Ptolemy (100-170 AD) explained these motions by introducing epicycles, or secondary circular motions. 


<img src="https://upload.wikimedia.org/wikipedia/commons/0/0b/PSM_V78_D326_Ptolemy.png">

However, this model was unable to accurately predict the motion of planets through time. 


<img src="https://physics.weber.edu/schroeder/ua/Epicycle.png">

<img src="https://study.com/cimages/multimages/16/cassini_apparent.jpg">


It took the work of a monk from the Baltic, Nicolaus Copernicus, to elaborate a much better model.
Copernicus realized that all the complexities of the system of Ptolemy could be eliminated by assuming that all planets revolve around the Sun. The next figure shows how this  new coperncian system explains the retrograde motion of  planets:

<img src="http://www.astronomy.ohio-state.edu/~pogge/Ast161/Unit3/Images/retrograde.gif">

However, even that was not good enough to track events such as ellipses or oppositions. The problem is that Copernicus still assumed perfect circular motion at constant speed.  

<a href="http://www.erasmatazz.com/library/book-reviews/book-reviews-for-2013/copernicus-and-modern-astro.html">This ruined the predictive power of his</a> otherwise excellent system.

At that point enters Kepler. He sits on a treasure: the most precise astronomical observations ever done, namely those of the danish astronomer Tycho Brahe, his master. With the methods explained here, Kepler will come up in his Astronomia Nova (1609) with a splendidly accurate reworking of the laws of astronomy. IN the next picture, you can see how shapes of the orbits have changed (from circular to elliptical, and how, contrary to Copernicus's asumption, the Sun is not the geometric center of the orbits. This drawing looks 'less symmetric' but is more accurate. The symmetries present in Nature do not always appear straightforward to the non-mathematicians.




<img src="http://www.erasmatazz.com/_Media/retromars-1_med_hr.jpeg">


With the material presented in this lecture, you can beat Ptolemy and Copernicus and accurately track any planet, for example Mars. This is indeed what Kepler endeavored to do. Here is his meticulous reconstruction of the trace of Mars, compared to a modern computer simulation:



These charts show the direction and distance of Mars relative to the Earth at the center, with oppositions and apparent retrograde motion approximately every 2 years and closest oppositions every 15–17 years due to Mars' eccentric orbit.


<img src="https://ars.els-cdn.com/content/image/1-s2.0-S1384107618303610-gr15.jpg">


Here is a recent reconstruction of the orbit of Mars seen from Earth. taken from wikipedia. Unlike the model of Ptolemy, you can see how the retrogradation loops vary in their shapes 

<img src="https://upload.wikimedia.org/wikipedia/commons/4/48/Mars_oppositions_2003-2018.png">






# Conclusion

I think it is fair to say that Kepler was the first major data scientist!
Like them, he gave an immense importance to data, the astronomic observations of Tycho Brahe. Like them, he spent much time - 8 years for the orbit of Mars- to unravel / mine the 'truth' contained in the data. He came out from this persistent labour with three fundamental laws of nature, something that data scientists usually do not succeed to do.
Like them he tried various theoretical models, of increasing accuracy, to best fit the data. Like them, he spent a lot of time doing numerical calculations, a great number thereof, and had to fill hundreds of pages of computations by hand as no computer was available at that time.




Kepler: the how but not the why

the first data scientist

The solid line shows the difference between Kepler's predictions and the observed position of Mars. (from Owen Gingerich, "Johannes Kepler and the Rudolphine Tables," Sky and Telescope, December, 1971, page 328)

http://ircamera.as.arizona.edu/NatSci102/NatSci102/images/marspos.gif

the why: Newton




file:///C:/Users/phili/OneDrive/Desktop/regressionlineaire/ACTIVITE-MARS/PlaneGeometry.pdf

```{r}

C   =c(0.1,-0.2)
S   =c(0,0)
gam =c(0,-1)
Line$new(A, B, extendA = TRUE, extendB = TRUE)
ell <- Ellipse$new(c(1,1), 5, 1, 30)
line <- Line$new(c(2,-2), c(0,4))
( Is <- intersectionEllipseLine(ell, line) )
ell$includes(Is$I1); ell$includes(Is$I2)

interddc <- function(S,gam,ang,C,R){
 
(x-C[1])^2+(y-C[2])^2  = R^2

SM vec gam = 0

  
  
  
return(c(x,y))
}


interdd(S,gam,ang,C,R)

```


plan de l'ecliptique

<iframe width="560" height="315" src="https://www.youtube.com/embed/KkiuYWuKTv0" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>


completely misleading

path of sun across the sky

<iframe width="560" height="315" src="https://www.youtube.com/embed/F95ArVB4Ye0" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>



the sun path along the ecliptic through the zodiac

<iframe width="560" height="315" src="https://www.youtube.com/embed/x1DqwC7lRpk" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>


<iframe width="560" height="315" src="https://www.youtube.com/embed/Xm13Kq_E1ik" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>


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



```{r}
getwd()
```

## Reconstruction of the orbit of the Earth

```{r}
mydata=read.csv("keplermars1.csv")
#summary(mydata)
head(mydata)

ls = pi *(mydata$lsdeg + mydata$lsmin/60)/180.
lm = pi *(mydata$lmdeg + mydata$lmmin/60)/180.

ls
lm

lmdeg0=334  
lmmin0=58
lsdeg0 =154 
lsmin0 =58

ls0=pi *(lsdeg0 + lsmin0/60)/180.
lm0=pi *(lmdeg0 + lmmin0/60)/180.
```


```{r}
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
yd
XY <- matrix(c(xd,yd),7,2,byrow=FALSE)
XY

plot(XY,xlim=c(-1,1),ylim=c(-1,1),asp=1)

sqrt(xd^2+yd^2)


```



```{r}
library(conicfit)
xd
yd

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



```



```{r}

XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=TRUE)
XY

XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=FALSE)
XY

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



plot(XY)

sqrt(xd^2+yd^2)



```





```{r}
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

XY <- matrix(c(1,7,2,6,5,8,7,7,9,5,3,7,6,2,8,4),8,2,byrow=FALSE)
XY

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



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


print('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
b
c
d
f



-xd^2+a*yd^2+b*xd*yd+c*xd+d*yd+f



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




```{r}
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


print('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
b
c
d
f


```


```{r}


xc=(2*a*c-b*d)/(4*a+b*b)
xc
yc=(b*c+2*d)/(-b*b-4*a)
yc
eps=f+xc*xc-a*yc*yc-b*xc*yc
eps


phi=atan(yc/xc)

# longitude of perigee

lonp=180.*(phi+lm0)/pi
lonp     # 283 IMCCE

# coord of perigee

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


print('delta')
delta=bb-4*aa*cc
delta


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

yc/xc


A=d2 # CP  demi grand axe 0.72
Ainv=1/d2 # IMCC 1.38


```



## Reconstruction of the orbit of the Mars



```{r}
mydata2=read.csv("keplermars2.csv")
#summary(mydata2)
mydata2

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} )$$



```{r}
xm=c()
ym=c()
```



```{r}

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


```{r}


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()
xm=c(xm,vv[1])
ym=c(ym,vv[2])
xm
ym




```




```{r}
xm
ym
```

```{r}
mydata2=read.csv("keplermars2.csv")
#summary(mydata2)
mydata2

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.



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))
}

```


Build the list of positions of Mars on its orbit:





```{r}
xm
ym
```


```{r}


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)
}



```


```{r}
  a=-0.9996
  b=-0.0002594
  c= 0.01495
  d=-0.01926
  f= 0.5242
gett1(1,1)
gett1(1,-1)


gett1(4,1)
gett1(4,-1)

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)
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)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')

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

vv2=gett2(1,-1)
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')








```

<table>
<tr>
<td><img height=250 src="imkk.png">
</td><td><img height=250 src="imkkk.png">
</td>
</tr>
</table>



```{r}

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))


vv1a=gett1(i,1)
lines(c(0,vv1a[1]),c(0,vv1a[2]),col='red')


vv1=gett1(i,-1)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')


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

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



vv2a=gett2(i,1)
lines(c(0,vv2a[1]),c(0,vv2a[2]),col='pink')

vv2=gett2(i,-1)
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

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

m1=vv2a
m2=m1+2*uu

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


```




```{r}

i=3
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))


vv1a=gett1(i,1)
lines(c(0,vv1a[1]),c(0,vv1a[2]),col='red')


vv1=gett1(i,-1)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')


t=lm1[i]-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')



vv2a=gett2(i,1)
lines(c(0,vv2a[1]),c(0,vv2a[2]),col='pink')

vv2=gett2(i,-1)
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

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

m1=vv2a
m2=m1+2*uu

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


```



```{r}

i=1
#
plot(xyDirect[,1],xyDirect[,2],type='l',
col='cyan',xlim=c(-1.5,1.5),ylim=c(-1.5,1.5),asp=1,main="b-b");par(new=TRUE)
lines(c(0,1),c(0,0))
vv1a=gett1(i,1)
lines(c(0,vv1a[1]),c(0,vv1a[2]),col='red')
vv1=gett1(i,-1)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')
vv2a=gett2(i,1)
lines(c(0,vv2a[1]),c(0,vv2a[2]),col='pink')
vv2=gett2(i,-1)
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

t=lm1[i]-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')

t=lm2[i]-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(-1.5,1.5),ylim=c(-1.5,1.5),asp=1,main="a-b");par(new=TRUE)
lines(c(0,1),c(0,0))
vv1a=gett1(i,1)
lines(c(0,vv1a[1]),c(0,vv1a[2]),col='red')
vv1=gett1(i,-1)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')
vv2a=gett2(i,1)
lines(c(0,vv2a[1]),c(0,vv2a[2]),col='pink')
vv2=gett2(i,-1)
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

t=lm1[i]-lm0
uu=c(cos(t),sin(t))
m1=vv1a
m2=m1+2*uu
lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')

t=lm2[i]-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(-1.5,1.5),ylim=c(-1.5,1.5),asp=1,main="b-a");par(new=TRUE)
lines(c(0,1),c(0,0))
vv1a=gett1(i,1)
lines(c(0,vv1a[1]),c(0,vv1a[2]),col='red')
vv1=gett1(i,-1)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')
vv2a=gett2(i,1)
lines(c(0,vv2a[1]),c(0,vv2a[2]),col='pink')
vv2=gett2(i,-1)
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

t=lm1[i]-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')

t=lm2[i]-lm0
uu=c(cos(t),sin(t))
m1=vv2a
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(-1.5,1.5),ylim=c(-1.5,1.5),asp=1,main="a-a");par(new=TRUE)
lines(c(0,1),c(0,0))
vv1a=gett1(i,1)
lines(c(0,vv1a[1]),c(0,vv1a[2]),col='red')
vv1=gett1(i,-1)
lines(c(0,vv1[1]),c(0,vv1[2]),col='green')
vv2a=gett2(i,1)
lines(c(0,vv2a[1]),c(0,vv2a[2]),col='pink')
vv2=gett2(i,-1)
lines(c(0,vv2[1]),c(0,vv2[2]),col='blue')

t=lm1[i]-lm0
uu=c(cos(t),sin(t))
m1=vv1a
m2=m1+2*uu
lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')

t=lm2[i]-lm0
uu=c(cos(t),sin(t))
m1=vv2a
m2=m1+2*uu
lines(c(m1[1],m2[1]),c(m1[2],m2[2]),col='black')




```



```{r}



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)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(1,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')

vv=gett1(2,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(2,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')


vv=gett1(3,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(3,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')


vv=gett1(4,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(4,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(5,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(5,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(6,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(6,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(7,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(7,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')



vv=gett1(8,1)
lines(c(0,vv[1]),c(0,vv[2]),col='red')
vv=gett1(8,-1)
lines(c(0,vv[1]),c(0,vv[2]),col='blue')


vv=gett2(1,1)
lines(c(0,vv[1]),c(0,vv[2]),col='pink')
vv=gett2(1,-1)
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))
```






```{r}
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)
}



```

```{r}


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
plot(xy)
```



```{r}


xym <- matrix(c(xm,ym),8,2,byrow=FALSE)
xym
plot(xym)
sqrt(xm^2+ym^2)

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


print('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
bm
cm
dm
fm



-xm^2+am*ym^2+bm*xm*ym+cm*xm+dm*ym+fm



```





```{r}
XY <- matrix(c(xm,ym),8,2,byrow=FALSE)
XY
plot(XY)
xy=XY
```



```{r}
XY <- matrix(c(xm,ym),8,2,byrow=FALSE)
XY
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


print('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
b
c
d
f

```

# The original work of Kepler

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


<img src="https://www.keplersdiscovery.com/Images/Keplers.Ch24%20Diagram.jpg">



date1,ls1degd,lm1degd,date2,ls2degd,lm2degd



```{r}
mydata2a=read.csv("keplermars2a.csv")
#summary(mydata2)
mydata2a

ls1 = pi *mydata2a$ls1degd/180.
ls2 = pi *mydata2a$ls2degd/180.
lm1 = pi *mydata2a$lm1degd/180.
lm2 = pi *mydata2a$lm2degd/180.


```



http://clea-astro.eu/lunap/Kepler/resolveuid/4183e5013da25f7eccf7b4a5abf345bf

http://clea-astro.eu/lunap/Kepler/resolveuid/9c1487b76a151b34e78f1b639d43bb59




Figure 1 : La construction de Kepler avec des paires de points.
Vous voyez sur la Figure 1 que par ces observations Kepler a pu donner la position de Mars pour chacune des paires de mesures, en supposant la Terre se déplaçant régulièrement sur une orbite circulaire (ce qui est une bonne approximation, mais Kepler avait pu aussi tracer la forme de l'orbite terrestre).
Avec la date on peut placer la Terre aux emplacements exacts sur son orbite autour du Soleil.

Regardez, par exemple, le point 5a. Il correspond au 10 mars 1585, une dizaine de jours avant l'équinoxe de printemps. Vu depuis 5a le Soleil est presque dans la direction du point g (le point gamma). La direction du Soleil vue depuis le point 2a (19 septembre) nous donne presque la direction du point g', opposé au point gamma.

En construisant les cinq positions de Mars pour les cinq paires de mesures, Kepler a vu que le Soleil n'était pas au centre de la trajectoire de Mars. L'excentricité (écart entre le centre réel et la position du Soleil, rapporté à la mesure du rayon maximal) est de 0.09, presque dix fois celle de la Terre (Si vous avez fait la construction, essayez de mesurer l'excentricité, vous serez surpris de la précision que l'on peut obtenir).
C'est ainsi que Kepler a pu faire l'hypothèse que la trajectoire de Mars était bien représentée par une ellipse.
Il a constaté aussi que la vitesse de déplacement des planètes n'était pas constante. Quand la planète est loin du Soleil sa vitesse est plus faible. Mais alors, me direz vous, c'est aussi vrai pour la Terre et le temps, mesuré par la position du Soleil, va être affecté d'une variation traduisant les changements réguliers de vitesse de la Terre sur son orbite.
Oui, c'est ce qu'on appelle l'équation du temps.




## 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$$


<img src="http://faculty.uca.edu/njaustin/PHYS1401/_Media/orbitMap01.jpg">


### 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

<a href="https://faculty.uca.edu/njaustin/PHYS1401/Laboratory/kepler.pdf">PHYS1401 Mars Lab</a>

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:

<img src="http://farside.ph.utexas.edu/Books/Syntaxis/Almagest/img617.png">


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

```{r}
paste('u = ',4)

mydata=read.csv("https://mathstat.dal.ca/~fullsack/DATA/pldata1.csv")
summary(mydata)
head(mydata)
```



## By measurements





The Earth's heliocentric ecliptic longitude is opposite to the Sun's geocentric ecliptic longitude. If the Earth's perihelion is at 103 degrees heliocentric longitude, then the geocentric longitude of the Sun at that time is 283 degrees, or 13 degrees past the December solstice.

comment Kepler a trob

file:///C:/Users/phili/OneDrive/Desktop/regressionlineaire/ACTIVITE-MARS/CLEA_CahiersClairaut_117_04.pdf


# LINKS

file:///C:/Users/phili/OneDrive/Desktop/regressionlineaire/VENUS/CLEA_CahiersClairaut_127_08.pdf

Les Cahiers Clairaut, bulletin du CLEA (comité de
liaison enseignants et astronomes) n° 117, printemps 2007, article de Blaise Simon n° 118, été
2007, Les orbites de la Terre et de Mars : la
première loi de Kepler, article de Béatrice Sandré.
KEPLER par Philippe Depondt et Guillemette de
Véricourt (Édition du Rouergue 2005).
KEPLER, astronome, astrologue par Gérard Simon
(Gallimard 1979).
http://acces.inrp.fr/clea/lunap/Kepler

https://faculty.uca.edu/njaustin/PHYS1401/Laboratory/kepler.html

<!--  art and mannerism at the Prague court


https://www.metmuseum.org/toah/hd/rupr/hd_rupr.htm

Described by a noted contemporary as “the greatest art patron in the world,” Rudolf II Habsburg (1552–1612), king of Hungary and Bohemia, and Holy Roman Emperor, raised court patronage in post-Renaissance Europe to a new level of breadth and extravagance. The thriving city and era over which he reigned, from 1583 until his death twenty-nine years later, is known as Rudolfine Prague. Seat of the emperor almost uninterruptedly from the mid-fourteenth century, Prague became, under Rudolf’s guidance, one of the leading centers of the arts and sciences on the continent. His taste for outstanding decoration and fantastic imagery were legendary, while his ambition and insight as a patron and collector changed the way art would be viewed by future generations.

Prague was founded in the late ninth century around the castle district of Hradcany, the oldest of the city’s four districts, and became the imperial residence during the reign of Charles IV in the 1300s. Rudolf II, son of Maximilian II, was named Holy Roman Emperor in 1576, and returned the court to Prague in 1583, after its temporary relocation to Vienna. As the city once again became the political and cultural focus of the empire and Hradcany earned the status of a royal town, Rudolf II brought into his service some of the most important European artists, architects, scientists, philosophers, and humanists, turning Prague into what has been referred to as a “Parnassus of the arts.” The astronomers Tycho Brahe (1546–1601), who was made Imperial Mathematician in 1599, and Johannes Kepler (1571–1630), who first served as Brahe’s assistant and then succeeded him in 1601, established observatories in Prague. The emperor’s ambitions as an architectural patron are evidenced by his redesign and expansion of the castle, the construction of a new town hall and archbishop’s palace, and the commissioning of several new churches. Yet, it was in painting, sculpture, and the decorative arts that Rudolf’s impact was most celebrated and distinct. Among the artists who came to the imperial court were the painters Bartholomeus Spranger (1997.93), Hans von Aachen, Pieter Stevens, and Roelandt Savery (1999.223); the miniaturists Joris Hoefnagel (63.200.4) and his son Jakob Hoefnagel; the sculptor Adriaen de Vries (41.190.534); the goldsmiths Paulus van Vianen and Wenzel Jamnitzer (17.190.620); and Aegidius Sadeler, who, as Imperial Printmaker (from 1597), popularized the emperor’s image (51.501.6469) and disseminated knowledge of his artists’ works. Foreign artists were especially prized by Rudolf because they gave international weight to his domain and satisfied his taste in art—for Italian and Netherlandish work, in particular, fostered at the Habsburg court in Spain, where he was educated.

Rudolf’s unrivaled passion for collecting culminated in one of the greatest of princely Kunstkammers, which contained small bronzes, works in cut stone, medallions and ivories, books and drawings, coins, scientific instruments and natural objects, as well as some paintings. Often dismissed as an unsystematic cabinet of curiosities intended for amusement or wonder, Rudolf’s Kunstkammer probably served as a place of refuge or serious contemplation, and indeed reflected, in an ordered way, the broader scientific and artistic interests of the court. Gerhard Emmoser’s dazzling celestial globe with clockwork (17.190.636), Female Nude (17.190.467), Apollo (41.190.534), Allegory of Virtues and Vices at the Court of Charles V (17.190.745), listed in the early seventeenth-century inventory of the collection, exemplify this close association. The intricate and colorful allegorical portraits of Rudolf II painted by Giuseppe Arcimboldo point to the direct collaboration that took place between artists and scholars. Rudolf also amassed a collection of paintings numbering in the thousands. Italian works were common, including those by Paolo Veronese (10.189), Correggio, and Leonardo da Vinci; as were works by northern European masters, including Albrecht Dürer and Pieter Bruegel the Elder.

In addition to stimulating culture in Prague through enlightened patronage and collecting, Rudolf II offered direct support of the arts, visiting his artists in their workshops and raising the status of the local painters’ guild from the level of a craft to that of a liberal art. Soon after the emperor’s death in 1612, his collections and court entourage were largely dispersed, leaving little in situ. The legacy of Rudolfine Prague, while apparent in some of the surviving buildings and landscaping of the old city, is best appreciated in the works of art commissioned and collected by one of Europe’s most influential and adventurous patrons.

-->

