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

1 files changed, 48 insertions, 9 deletions

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)