# Aurora Hiveley, 04/11/24, Assignment 22

Help:= proc(): 
print(` DO NOT MODIFY PARAMETERS, JUST READ FILE IN MAPLE WITH SEMICOLON`): 
print(` MAKE SURE TO CLICK PLAY ON SECOND PLOT TO VIEW ANIMATION`):
end:


### problem 2
# F_G = GMm/d^2 where d is the distance between earth and some other heavenly body
# then d/dr(F_g) = -2GMm/d^3, so the change in gravitational force with distance 
# (responsible for the tides as the difference between F_G on the "front" and "back"
# of the planet is what causes tides) goes as 1/d^3. 
# since the lunar tides are known to be stronger than the solar tides, we conclude that dF_G/dr
# from the moon is stronger than from the sun, and so m/d^3 is larger for the moon
# than for the sun. since the moon and the sun appear to be the same size from earth,
# we might conclude that the moon and sun are the same distance d from earth.
# then the discrepancy in dF_G/dr can only be explained by the moon being more massive
# than the sun. since the two appear to have the same size (volume), we have that 
# mass_moon > mass_sun implies density_moon > density_sun (as density = mass/volume).



### problem 3: project checkpoint
# pablo wrote some great code which generates the transition matrix for all english words 
# of length >=k. using this transition matrix code, i generated some words which have a 
# high likelihood of being english words, but actually aren't.

# my code generateWord(K,A) takes in a word length K and a starting letter A (input should be
# lowercase letter, no quotes or parens) and, for each successive letter of the word, selects 
# the letter which most often follows the previous letter. in other words, when the current letter 
# is "i", we identify the maximum entry in row "i". say this is in column "j", then letter "j" is 
# next in the word (after adjusting indices for pablo's "start" and "end" rows).
# lastly, we check whether or not the words are in ENG(), and returns them if not


## some example calls and results
# generateWord(5,c); -- > [c, o, n, g, e]
# generateWord(7,o); --> [o, n, g, e, r, e, r]
# generateWord(4,q); --> [q, u, n, g]



### problem 4: eclipse diagram
with(plots): 
with(plottools):


## some helper procs
# calculate slope and intersect of line through two points
lineareq := proc(p1,p2) local s,b:
s := (p2[2]-p1[2])/(p2[1]-p1[1]):
b := p1[2] - p1[1]*s:
[s,b]:
end:

# find point of intersection of two lines determined by two points each
# one line through p1 and p2, the other line through p3 and p4
poi := proc(p1,p2,p3,p4)  local x,y,sa,sb,ba,bb:
sa := lineareq(p1,p2)[1]:
ba := lineareq(p1,p2)[2]:
sb := lineareq(p3,p4)[1]:
bb := lineareq(p3,p4)[2]:

x := fsolve(sa*z+ba = sb*z+bb, z=0..10):
y := sa*x + ba:
[x,y]:
end:


## plot solar eclipse diagram
SE := proc(scent, srad, mcent, mrad, ecent, erad): 
local S,M,E,Eo,Mo,a,sp1,sp2,mp1,mp2,l1,l2,l3,l4,l5,l6,L,P,P1,P2,shadow:

# shapes for sun, moon, and earth
S := plot([scent[1]+srad*cos(t), scent[2]+srad*sin(t), t = 0..2*Pi], axes=none,filled=true, color=yellow):
M := plot([mcent[1]+mrad*cos(t), mcent[2]+mrad*sin(t), t = 0..2*Pi], axes=none,filled=true,color=gray):
E := plot([ecent[1]+erad*cos(t), ecent[2]+erad*sin(t), t = 0..2*Pi], axes=none,filled=true,color=blue):

# orbits of earth and moon
Eo := circle(scent,sqrt(ecent[1]^2 + ecent[2]^2),color=lightblue,linestyle="dash"):
Mo := circle(ecent,sqrt((mcent[1]-ecent[1])^2 + (mcent[2]-ecent[2])^2),color=lightgrey,linestyle="dash"):

## we do some geometry
# angle between sun and moon
a := arctan((mcent[2]-scent[2])/(mcent[1]-scent[1])):
# equatorial points on sun which emit rays
sp1 := scent + [srad*sin(a), -srad*cos(a)]:
sp2 := scent + [-srad*sin(a), srad*cos(a)]:
# equatorial points on moon which intersect sun rays
mp1 := mcent + [mrad*sin(a), -mrad*cos(a)]:
mp2 := mcent + [-mrad*sin(a), mrad*cos(a)]:

# sun rays between sun and moon equators
l1 := line(sp1,mp1,color=yellow):
l2 := line(sp1,mp2,color=yellow):
l3 := line(sp2,mp1,color=yellow):
l4 := line(sp2,mp2,color=yellow):

# point of intersection of the two more angled rays
P := poi(sp1,mp1,sp2,mp2): # umbra
# shadow rays between moons equator and umbra
l5 := line(mp1,P,color=grey):
l6 := line(mp2,P,color=grey):

# two points which estimate the size of the penumbra 
# note that only one is *truly* necessary, but two ensure that we did it correctly
P1 := poi(sp1,mp2,ecent,[0,ecent[2]]):
P2 := poi(sp2,mp1,ecent,[ecent[1]-0.001,0]):

# penumbra
shadow := disk(P,sqrt((P[1]-P1[1])^2 + (P[2]-P1[2])^2), color=black):

# put everything on one plot
L := [l1,l2,l3,l4,l5,l6]:
display(S,M,E,Eo,Mo,shadow,pointplot(P),L,'view' = [-2..1.38*ecent[1],-2..1.38*ecent[2]]);
end:


# sample call
SE([0,0],1.2,[7,7],0.16,[8.2,8.2],0.8);



### problem 5: challenge animation


# NOTE: FOR WHATEVER REASON, REWRITING THIS AS A PROC DOESN'T DISPLAY THE WHOLE PLOT 
# WHEN CODE IS RUN. SO THIS IS NOT A PROC, BUT THE INPUT PARAMETERS BELOW ARE CHANGEABLE, 
# AND READING THIS .TXT FILE SHOULD CREATE PLOT.

# DON'T FORGET TO CLICK PLOT AND PRESS PLAY TO SEE THE ANIMATION

earthPeriod := 365.5:
moonPeriod:= 27.32: 

with(plots):
with(plottools):

# rescale periods
# earthPeriod = 2*Pi
# earthPeriod/2*Pi = moonPeriod/per*(2*Pi)
per := moonPeriod/earthPeriod:

sun := pointplot([[0,0]],color=yellow,symbol=solidcircle,symbolsize=80):

# position of earth at time t for t from 0 to 2pi
ep := proc(t)
[0.5*cos(t), 0.5*sin(t)]:
end:

# position of moon at time t for t from 0 to 2pi
mp := proc(p,t)
[p[1] + 0.2*cos(t), p[2]+0.2*sin(t)]:
end:

# define earth plot at position (x,y)
earth := proc(x,y) 
pointplot([[x,y]],color=blue,symbol=solidcircle,symbolsize=40):
end proc:

# define moon plot at position (x,y)
moon := proc(x,y) 
pointplot([[x,y]],color=grey,symbol=solidcircle,symbolsize=20):
end proc:


# animate orbits of earth and moon
ae := animate(earth,ep(t),t=0..2*Pi, scaling=constrained, axes = none, frames = 100):
am := animate(moon,mp(ep(t),t/per),t=0..2*Pi,scaling=constrained, axes = none, frames = 100):

# assumed circular orbit of earth
eo := circle([0,0], 0.5, color=lightblue):

# display all items together
display(sun,eo,ae,am);



### a note for dr z: to animate rotation of earth, i recommend using the rotate() command
## something like animate(display, rotate(earth, t), t=0..2*Pi):
## this won't work on pointplot (maple doesn't care about rotating points)
## but if you encoded the earth another way, this should work