From 61a5bbdaf45bf29f8da9e7657bb5b23a575ef27c Mon Sep 17 00:00:00 2001 From: Bernhard Tittelbach Date: Wed, 21 Mar 2012 23:56:04 +0000 Subject: added moon, added example constellations, added timestep choosen for shannon limit git-svn-id: https://svn.spreadspace.org/mur.sat@327 7de4ea59-55d0-425e-a1af-a3118ea81d4c --- tools/photodiodesim/pd_sim.py | 57 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 48 insertions(+), 9 deletions(-) (limited to 'tools') diff --git a/tools/photodiodesim/pd_sim.py b/tools/photodiodesim/pd_sim.py index 4416971..14cbd22 100755 --- a/tools/photodiodesim/pd_sim.py +++ b/tools/photodiodesim/pd_sim.py @@ -17,13 +17,31 @@ w_y = np.random.ranf() * np.pi/2 #rad/s w_z = np.random.ranf() * np.pi/2 #rad/s filename_export_mat = os.path.join(os.path.curdir,"pd_sim.mat") filename_export_csv = os.path.join(os.path.curdir,"pd_sim.csv") +## fun constellations to simulate: +# sun and earth close next to each other and no light from moon sun = [1,1,1] +earth = [1,8,-5] +moon = [1,1,1] +# sun and earth next to each other, moon at maximum +#~ sun = [1,1,1] +#~ earth = [1,8,-5] +#~ moon = [-1,-1,-1] +# sun hidden behind earth, only moon visible +#~ sun = [1,1,1] +#~ earth = [1,1,1] +#~ moon = [1,-1,0] +#satellite between sun and earth +#~ sun = [1,1,1] +#~ earth = [-1,-1,-1] +#~ moon = [1,-1,0] +#moon and earth opposite to sun +#~ sun = [1,1,1] +#~ earth = [-1,-2,-1] +#~ moon = [-1,2,-1] sun_intensity = 1.0 -earth = [0,-1,-1] earth_reflection_intensity = 0.10 * sun_intensity # 0 to turn earth off ;-) moon_reflection_intensity = 0.07 * sun_intensity # 0 to turn moon off ;-) orbit_height = 310e3 #m -numsamples = 10000 numperiods = 4 sensors = np.matrix([ [ 1,0,0], @@ -37,8 +55,9 @@ sensors = np.matrix([ #calculate simduration so we see numperiods periods of the smalest angular speed w_xyz = np.matrix([w_x,w_y,w_z]) duration_full_period = np.prod(np.pi / np.array(filter(lambda x: x>0.0,w_xyz))) if w_x or w_y or w_z else 1.0 +minimum_period = np.min(np.pi / np.array(filter(lambda x: x>0.0,w_xyz))) if w_x or w_y or w_z else 1.0 simduration = numperiods* duration_full_period #in seconds -simtimestep = simduration / numsamples #s +simtimestep = minimum_period / 10.0 #s #calculate angle between normalvector of earthdisc and edge of earthdisc as seen from satellite earth_radius = 6e6 #m earth_distance = earth_radius + orbit_height #m @@ -48,19 +67,25 @@ simtimeline = np.arange(0,simduration,simtimestep) angular_speed_sum_vectors = ( w_xyz.transpose() * simtimeline ).transpose() % (2 *np.pi) sun = sun / np.linalg.norm(sun) earth = earth / np.linalg.norm(earth) +moon = moon / np.linalg.norm(moon) angle_sun_earth = np.arccos(max(min(np.dot(sun,earth),1.0),-1.0)) +angle_sun_moon = np.arccos(max(min(np.dot(sun,moon),1.0),-1.0)) +angle_earth_moon = np.arccos(max(min(np.dot(earth,moon),1.0),-1.0)) +print angle_sun_earth, angle_sun_moon, angle_earth_moon #print angle_sun_earth /np.pi*180, angle_earthdisc_visible /np.pi*180 print "Simulation Input:" -print "　 Duration:", simduration -print "　 Step Size:", simtimestep +print "　 Duration : %.4f [s]" % simduration +print "　 Step Size: %.4f [s]" % simtimestep print "　 Angular Speed x-axis: %.3f [rad/s] period: %.3f [s]" % (w_x, np.pi / w_x) print "　 Angular Speed y-axis: %.3f [rad/s] period: %.3f [s]" % (w_y, np.pi/ w_y) print "　 Angular Speed z-axis: %.3f [rad/s] period: %.3f [s]" % (w_z, np.pi/w_z) print "　 Full period %.3f [s]" % duration_full_period print " Earth " + ("Light intensity: %f" % earth_reflection_intensity if earth_reflection_intensity else "is switched OFF") +print " Moon " + ("Light intensity: %f" % moon_reflection_intensity if moon_reflection_intensity else "is switched OFF") print " Sphere surface, covered by Earth as seen from satellite: %.2f %%" % ((earth_radius**2) / (4.0 * earth_distance**2) * 100.0) -print " Satellite " + (u"Is (!) " if angle_sun_earth < angle_earthdisc_visible else u"is not") + u" in Earth's shadow (Angle Sun-Earth: %d°)" % (angle_sun_earth /np.pi*180.0) +print " Satellite " + (u"Is (!)" if angle_sun_earth < angle_earthdisc_visible else u"is not") + u" in Earth's shadow (Angle Sun-Earth: %d°)" % (angle_sun_earth /np.pi*180.0) +print " Moon " + (u"Is (!)" if angle_earth_moon < angle_earthdisc_visible else u"is not") + u" hidden behind the Earth (Angle Earth-Moon: %d°)" % (angle_earth_moon /np.pi*180.0) print "" @@ -114,6 +139,20 @@ def earthdisc_light_rel_intensity(angle): earth_model = earthdisc_light_rel_intensity +def moon_reflection_from_sun_rel_intensity(angle): + """ reflected light from moon is a function of angle between sun and moon as seen from satellite if < 90% + also: moon can be hidden by earthdisc + """ + global angle_sun_moon, sensor_model, angle_earth_moon + if angle_earth_moon < angle_earthdisc_visible: + return 0.0 + elif angle_sun_moon < np.pi/2.0: + return angle_sun_moon / (np.pi/2.0) * sensor_model(angle) + else: + return sensor_model(angle) + +moon_model = moon_reflection_from_sun_rel_intensity + def _quote(str): return '"'+ re.sub(r'(?!\\)"','\\"',str) +'"' @@ -129,13 +168,13 @@ suntrace=[] earthtrace=[] for cur_rot in angular_speed_sum_vectors: rot_sensors=rotateVector(sensors,cur_rot) - #sinphi=np.cross(sens1,sun) # /1 / 1 angles_to_sun=np.array(np.arccos(np.dot(rot_sensors,np.matrix(sun).T))).flatten() angles_to_earth=np.array(np.arccos(np.dot(rot_sensors,np.matrix(earth).T))).flatten() - #print [sensor_rel_sensitivity_cos(a) for a in angle_sun] + angles_to_moon=np.array(np.arccos(np.dot(rot_sensors,np.matrix(moon).T))).flatten() light_intensity_per_sensor = \ sun_intensity * np.array([sun_model(a) for a in angles_to_sun]) + \ - earth_reflection_intensity * np.array([earth_model(a) for a in angles_to_earth]) + earth_reflection_intensity * np.array([earth_model(a) for a in angles_to_earth]) + \ + moon_reflection_intensity * np.array([moon_model(a) for a in angles_to_moon]) senstrace.append(light_intensity_per_sensor) suntrace.append(angles_to_sun) -- cgit v1.2.3