热度 3||
# -*- coding: utf-8 -*-
"""
Created on Fri Jul 8 21:00:50 2018
@author: YDX
"""
from sympy import symbols, cos, sin, acos, sqrt, pprint,diff,solve,nsolve, Eq, latex, integrate,pi, simplify, Matrix
from sympy.plotting import plot
import matplotlib.pyplot as plt
import numpy as np
import sympy as sym
# alpha -- angle between sun ray and tilt axis
# tilt -- titl axis
# phis -- angle around sun
tilt, phis, d, lat, sp, t = symbols(' tilt phis d lat sp t')
# z^ = cos(tilt) z'^ - sin(tilt) y'^
# y^ = sin(tilt) z'^ + cos(tilt) y'^
# x^ = x'^
# sun ray: [sin(phis) x'^, - cos(phis) y'^,0]
#point: [cos(lat)*cos(d), cos(lat)*sin(d), sin(lat)]
m = sym.zeros(3,3)
m[0,0]=1
m[1,2] = sin(tilt)
m[1,1] = cos(tilt)
m[2,1]= -sin(tilt)
m[2,2]= cos(tilt)
print(m)
loc = Matrix([cos(lat)*sin(d), -cos(lat)*cos(d), sin(lat)])
east = -loc.cross(Matrix([0,0,1]))
tnorth=-Matrix([-sin(lat)*sin(d), sin(lat)*cos(d), cos(lat)])
print ("east", east)
east_len = simplify(sqrt((sym.Transpose(east)*east)[0,0]))
print (east_len)
east = east/east_len
print ("east norm", east)
east2 = simplify( loc.cross(tnorth))
print ("east 2", east2)
north = Matrix([0,0,1])
ray = Matrix([sin(phis), - cos(phis),0])
loc_sun = sym.Transpose(m) * loc
east_sun = sym.Transpose(m) * east2
tnorth_sun = sym.Transpose(m) * tnorth
dotp = sym.Transpose(ray)*loc_sun
dotn = sym.Transpose(ray)*tnorth_sun
dote = sym.Transpose(ray)*east_sun
#print(dotp.shape, dotp)
dot=simplify(dotp[0,0])
print("dotn=",dotn)
dotn=simplify(dotn[0,0])
print("dotn=",dotn)
sp = acos(dot)
print("angle", sp)
sp2= acos(dotn)
print(latex(sp))
dlen= 23+56/60+4/3600
ylen = 365.2422
oe = 2*np.pi * 24/dlen
os = 2*np.pi/ylen
d2rad = np.pi/180
til = 23.43688*d2rad
lati = 40*d2rad
def angle(t, tilt, lat):
d = oe *t
phis = os *t
return np.arccos((np.cos(d)*np.cos(lat)*np.cos(tilt) + np.sin(lat)*np.sin(tilt))*np.cos(phis) + np.sin(phis)*np.sin(d)*np.cos(lat))
ts=np.arange(0,365,0.005)
angles = angle(ts, til, lati)
plt.plot(ts[angles>np.pi/2], angles[angles>np.pi/2])
spd =sp.subs({lat:lati, tilt:til, phis:os*t, d:oe*t})
sp2d =sp2.subs({lat:lati, tilt:til, phis:os*t, d:oe*t})
print(latex(spd))
sols=[]
morning_angles=[];
day=0
lasts=0
intv = 0.2
a1 = np.arccos(-np.tan(til)*np.tan(lati))
d1=nsolve(Eq(spd,np.pi/2), t, 0)
d2=nsolve(Eq(spd,np.pi/2), t, d1+0.4)
d1=nsolve(Eq(spd,np.pi/2), t, 0)
d3=nsolve(Eq(spd,np.pi/2), t, 1)
print (d1,d2,d3, d1*24, d2*24, (d2-d1)*24, (d3-d1)*24)
d5= d1*oe
print ("diff=", a1, a1/d2rad, d5, (d5-a1)*24*60/oe)
lasts=0
#lasts = nsolve(Eq(spd,np.pi/2), t, [lasts+0.01,lasts+intv], solver='bisect')
#print (lasts)
#raise SystemExit
while lasts<185:
try:
lasts = nsolve(Eq(spd,np.pi/2), t, [lasts+0.0001,lasts+intv], solver='bisect')
sols.append(lasts)
ma = sp2d.subs({t:lasts}) /d2rad
sols.append(ma)
intv = 0.4
#print ("Found " +str(lasts))
except:
#print ("no solution", lasts, lasts+intv)
lasts = lasts+intv
days = np.array(sols)
sz = days.size
if(sz%4):
sz = sz- sz%4
nd=np.reshape(days[0:sz], (-1,4))
print(nd)
Powered by Discuz! X2.5