# Modeling dropping an object from a height

g = 2 #acceleration of "gravity"
y = 0 #starting y value.  As we move downward, y will INCREASE.  The variable y is POSITION
v = 0 #the velocity starts at zero, since we are just "letting go" of the object
t = 0 #starting time
dt = 0.1 # this is our "delta t" which should be SMALL as we approximate the object's motion and take "snapshots" every dt

# we are going to store the points (t, y) in a list to plot later
data = [(t,y)]  #so the data holds just [(0,0)] to start.

#Now we make a table of t, y, v.
nsteps = 200 #the number of snapshots.  So we will go for nsteps*dt time units.
print("t      y      v")  # the top row

for k in range(nsteps):    #this means that k goes from 0 until nsteps-1 and we loop
    #print n(t, digits =3), n(y, digits =3), n(v, digits =3)   #self explanatory (digits tells it how many decimal places to display)
    # now we iterate for the next snapshot
    t = t + dt            # time advances
    y = y + v*dt          # the "delta y" is equal to v*dt, since distance = rate times time!
    v = v + g*dt          # same idea
    data.append((t,y))    # add the new point to the list of points
list_plot(data)           # when the table is done, plot the t vs. y graph

# Modeling population growth
#
#  we will assume that P'(t) = aP(t), where a is a smallish constant

a = 0.2 #constant
P = 100 #starting P value.   
v = P*a #the starting velocity in people per time units (say, years)
t = 0 #starting time
dt = 0.1 # this is our "delta t" which should be SMALL as we approximate the object's motion and take "snapshots" every dt

# we are going to store the points (t, y) in a list to plot later
data = [(t,P)]  #so the data holds just [(0,0)] to start.

#Now we make a table of t, y, v.
nsteps = 200 #the number of snapshots.  So we will go for nsteps*dt time units.

print("t      P      v")  # the top row

for k in range(nsteps):    #this means that k goes from 0 until nsteps-1 and we loop
    #print n(t, digits =3), n(P, digits =3), n(v, digits =3)   #self explanatory (digits tells it how many decimal places to display)
    # now we iterate for the next snapshot
    t = t + dt            # time advances
    P = P + v*dt          # the "delta P" is equal to v*dt, since distance = rate times time!
    v = P*a               # same idea
    data.append((t,P))    # add the new point to the list of points
list_plot(data)           # when the table is done, plot the t vs. y graph

#catenary simulation
# we will try to create an approximate solution to the differential equation
# tan(t) = s, where t is theta and s is "ell" (this is a standard letter, and "ell" is hard to read in this font.)

t=0   #start at the bottom of the chain
x=0
y=0
ds =0.01 # this is what we increment length by; a smaller value yields a better approx
s = 0 #the starting arc length from (0,0).  It will increase as we move to the right.
nsteps = 300  #self explanatory

data=[(x,y)]  # So we start with (0,0) and then we will append and finally plot
for k in range(nsteps):   #the expression "range(nsteps)" just means, "the numbers from 0 to nsteps-1, inclusive"
    #we will move by a length of ds
    # in the direction with slope = tan(t)
    # and tan(t)=s, the total distance so far
    x += ds*cos(t)     #another programming expression: "x += u" means "increase x by the value u;" i.e. x = x+u
    y += ds*sin(t)
    
    #include the printout for debugging
    #print k, x, y,t
    #append the point
    data.append((x,y))
    
    #now increment
    s += ds    #computing the total length so far
    t = atan(s)  # the angle t is such that tan(t) = s; i.e., the SLOPE equals s at this point
    
#loop is over, now plot
list_plot(data,aspect_ratio=1)

# compare this with a "real" catenary
  #you need to state that x is a variable used in the plot command

%var x
approx = list_plot(data)
real = plot((exp(x)+exp(-x))/2-1, (x,0,1.8), color='red')    #exp(x) means e^x
show(approx+real,aspect_ratio=1)

#square-wheeled trike simulation
# we will try to create an approximate solution to the differential equation
# 1- tan(t) = s, where t is theta and s is "ell" (this is a standard letter, and "ell" is hard to read in this font.)

t=pi/4   #start at the left of the surface, with an angle of 45 degrees (square is on its corner)
x=-.8814 #this is the x-coordinate where the actual curve starts.  We could start with other values if we wanted.
y=0
ds =0.1 # this is what we increment length by; a smaller value yields a better approx
s = 0 #the starting arc length from (0,0).  It will increase as we move to the right.
nsteps = 30  #self explanatory

data=[(x,y)]  # So we start with (0,0) and then we will append and finally plot
for k in range(nsteps):   #the expression "range(nsteps)" just means, "the numbers from 0 to nsteps-1, inclusive"
    #we will move by a length of ds
    # in the direction with slope = tan(t)
    # and tan(t)=s, the total distance so far
    x += ds*cos(t)     #another programming expression: "x += u" means "increase x by the value u;" i.e. x = x+u
    y += ds*sin(t)
    
    #include the printout for debugging
    #print k, x, y,t
    #append the point
    data.append((x,y))
    
    #now increment
    s += ds    #computing the total length so far
    t = atan(1-s)  # the angle t is such that 1-tan(t) = s; i.e., the SLOPE equals s at this point
    
#loop is over, now plot
list_plot(data,aspect_ratio=1)

def spiralOfcirc(phase, twist, ratio, rmin, ncircles):
    r=rmin
    g=1+ratio
    sc=Graphics()
    xc=0; yc=0; theta=phase;radius=0
    for _ in range(ncircles):
        sc += circle((xc,yc),r,fill=True, facecolor='white')
        
        theta +=twist
        
        radius += r*g
        xc += g*r*cos(theta); yc +=g*r*sin(theta)
        r *= ratio
        
    return sc

show(spiralOfcirc(0,pi/6,1.1,1,12))

r=1.01
g=(1+sqrt(5))/2
ph=2*pi/g
flower=[spiralOfcirc(t*ph,pi/12,r^(t+1),1,12) for t in range(12)]
spirals=Graphics()
for s in flower:
    spirals += s
spirals.show(axes=false,xmin=-30,xmax=30,ymin=-30,ymax=30)
save(spirals,"spirals.png",axes=false,xmin=-30,xmax=30,ymin=-30,ymax=30)

sp=spirals(axes=false)

spirals=Graphics()
n=2
mult=1.1
rmin=.5
ncircs=12
phase=0
turn = 2*pi/12
for k in range(n):
    
    spirals += spiralOfcirc(phase,turn,mult,rmin,ncircs )
    phase += 2*turn
    #rmin *= mult
spirals.show()

spirs=Graphics()
for f in frames2:
    spirs +=f
show(spirs)

nSpirs=3
nTwists=6
startRadius=1
ratio=1.05
multiplier=1/ratio

spiralAn=[]
size=5
for tw in range(nTwists):
    r=startRadius
    spiral=Graphics()
    #spiral += polygon([(-size,-size), (size,-size), (size,size),(-size,size)])
    spiral += circle((0,0),size, fill=False)
    for t in range(nSpirs-1):
    
        spiral += spiralOfcirc(tw*2*pi/nTwists+t*2*pi/nSpirs ,pi/12,ratio,r,nSpirs)
        r *= multiplier
    spiralAn.append(spiral)

frames2=[spiralOfcirc(t*pi/6 ,pi/12,1.2,1,10) for t in range(12)]
spirs=Graphics()
for f in frames2:
    spirs +=f
show(spirs)

frames2=[spiralOfcirc(t*2*pi*0.388 ,pi/12,1.01^t,1,10) for t in range(15)]
spirs=Graphics()
for f in frames2:
    spirs +=f
show(spirs)

frames2=[spiralOfcirc(t*2*pi*0.388 ,pi/12,1.1,1.05^t,10) for t in range(18)]
spirs=Graphics()
for f in frames2:
    spirs +=f
show(spirs)

ph=pi*(1+sqrt(5))
n= 100
radius = .5
mult = 1.00
pic=Graphics()
for t in range(n):
    r=sqrt(t)
    a = t*ph
    #radius = .4+sqrt(t)*.025
    #pic += point((r*cos(a),r*sin(a)))
    pic += circle((r*cos(a),r*sin(a)),radius,fill=True, facecolor=hue(t/n,.2,.9))
max=20   
show(pic,aspect_ratio=1,axes=false,xmin=-max, xmax=max, ymin=-max,ymax=max)
save(pic,"spirals-g.png",aspect_ratio=1,axes=false,xmin=-max, xmax=max, ymin=-max,ymax=max)

show(pic,aspect_ratio=1)

%var r, t
n=13
f=2*pi/13
sp1=Graphics()
for k in range(n):
    u=f*k
    sp1 += polar_plot(e^(t+u),(t,-u,pi-u),aspect_ratio=1)

n=8
f=2*pi/8
sp2=Graphics()
for k in range(n):
    u=f*k
    sp2 += polar_plot(e^(pi-t+u),(t,u,pi+u),aspect_ratio=1,color='red')

    
show(sp2+sp1)

show(sp2+sp1)
save(sp1+sp2,"spir-fib.png",aspect_ratio=1, axes=false)

#transformations needed to build Edmark's tiles



def rot(t,vec):
    #rotate vector by angle t counterclockwise
    c = cos(t)
    s = sin(t)
    r= matrix([[c,-s],[s,c]])
    return r*vec

def trans(U, vec):
    #trans vec by vector U
    return vec+U

def cloneDown(A,B,C,D):
    #starting with quad ABCD with AB on "bottom", create quad A'B'C'D' below where D'=A, C'=B
    h=vector((1,0))
    DC = C-D
    AB = B-A
    r = AB.norm()/DC.norm()  #dilation ratio
    U=r*(A-D)
    #dilate, centered at A
    A1= A
    B1= r*(B-A) +A
    C1 = r*(C-A)+A
    D1=r*(D-A)+A
    # now translate so D1 moves to A
    A2=trans(U,A1)
    B2=trans(U,B1)
    C2=trans(U,C1)
    D2=trans(U,D1)
    #now rotate by t = angle to move CD to be parallel with AB
    #we will need a sign of a determinant to know direction
    m=Matrix([DC,AB])
    t=sgn(m.det())*acos(DC.dot_product(AB)/(DC.norm()*AB.norm()))
    A3=rot(t, A2-A)+A
    B3=rot(t, B2-A)+A
    C3=rot(t, C2-A)+A
    D3=rot(t, D2-A)+A
    
    return A3,B3,C3,D3

def cloneLeft(A,B,C,D):
    #starting with quad ABCD with AB on "bottom", create quad A'B'C'D' on left where B'=A, C'=D, A'D' on left
    X, Y, Z, W = cloneDown(D, A, B, C)
    return Y, Z, W, X

def cloneRight(A,B,C,D):
    #starting with quad ABCD with AB on "bottom", create quad A'B'C'D' on right where A'=B, D'=C, B'C' on right
    X, Y, Z, W = cloneDown(B,C,D,A)
    return W,X,Y,Z

#test run, using the coord of the "mother tile"
A=vector((0,0))
B=vector((1,0))
C= vector((1.1056158046757, 1.3611084974398))
D =vector((-0.294552318493650,1.0533128073132))

pic = Graphics()

pic += line([A,B,C,D,A])
X,Y,Z,W = cloneDown(A,B,C,D)
pic += line([X,Y,Z,W,X])
X,Y,Z,W = cloneLeft(A,B,C,D)
pic += line([X,Y,Z,W,X])

X,Y,Z,W = cloneRight(A,B,C,D)
pic += line([X,Y,Z,W,X])
     
show(pic,aspect_ratio=1)

#draw full tiling

A=vector((0,0))
B=vector((1,0))
C= vector((1.1056158046757, 1.3611084974398))
D =vector((-0.294552318493650,1.0533128073132))

#the four coords below "almost work"
#A=vector((0,0))
#B=vector((1,0))
#C= vector((1, 1.5))
#D =vector((-.5,1))


nsteps =21 #how many levels

pic = Graphics()
pic += line([A,B,C,D,A])
pic += polygon([A,B,C,D],color=hue(1)) #mother quad (so filled in)
quadsRight=[(A,B,C,D)]
for row in range(nsteps):
    
    X,Y,Z,W = quadsRight[row]
    X,Y,Z,W = cloneDown(X,Y,Z,W)
    #pic += line([X,Y,Z,W,X])
    pic += polygon([X,Y,Z,W], color=hue(.5*int(mod(row+col+1,2))))
    quadsRight.append((X,Y,Z,W))
for row in range(nsteps):
    X,Y,Z,W = quadsRight[row]
    for col in range(1, row+3):
        X,Y,Z,W = cloneLeft(X,Y,Z,W)
        #pic += line([X,Y,Z,W,X])
        pic += polygon([X,Y,Z,W], color=hue(.5*int(mod(row+col,2))))

for row in range(nsteps):
    X,Y,Z,W = quadsRight[row]
    for col in range(1, 4):
        X,Y,Z,W = cloneRight(X,Y,Z,W)
        #pic += line([X,Y,Z,W,X])
        pic += polygon([X,Y,Z,W], color=hue(.5*int(mod(row+col,2))))



show(pic,aspect_ratio=1,axes=false)  
save(pic,"edmark.pdf",aspect_ratio=1,axes=false)  

row

.5*int(mod(5,2))

a=5
n=12
data =[]
for k in range(n):
    print k, a
    data.append([k,a])
    a = a^2
    
latex(table(data))

a=5
n=8
data =[]
for k in range(n):
    print k, a^2-a
    data.append([k,a^2-a])
    a = a^2
    
latex(table(data))

# find x,y such that x^2+y^2=p for primes p
def sumOfSquares(n):
    # find x,y s.t. x^2+y^=n or if fail, print "no"
    output = "no"
    u=floor(sqrt(n))
    for x in range(u+1):
        y= n-x^2
        if y.is_square():
            output = (x,sqrt(y))
    return output
#start value

p=2

#end
biggest = 1000
while p<biggest:
    print p, sumOfSquares(p)
    
    #iterate to next prime
    p = p.next_prime()
    

def sumOfSquares3(n):
    # find x,y,z such that x^2+y^2+z^2= n  or if fail, print "no"
    output = "no"
    u=floor(sqrt(n))
    for x in range(u+1):
        d = n-x^2
        if sumOfSquares(d) != "no":
            y,z = sumOfSquares(d)
            output = x,y,z
    return output

for k in [4..100]:
    print k, sumOfSquares3(k)

%var x,y
plot3d((x-y)^2+(sqrt(2-x^2)-9/y)^2,(x,0,sqrt(2)),(y,2,3))