added earthdisc to photodiode simulator

git-svn-id: https://svn.spreadspace.org/mur.sat@324 7de4ea59-55d0-425e-a1af-a3118ea81d4c

1 files changed, 77 insertions, 34 deletions

diff --git a/tools/photodiodesim/pd_sim.py b/tools/photodiodesim/pd_sim.py
index 0d95522..66cdacc 100755
--- a/tools/photodiodesim/pd_sim.py
+++ b/tools/photodiodesim/pd_sim.py
@@ -17,37 +17,51 @@ w_z = np.random.ranf() * np.pi/2 #rad/s
 filename_export_mat = "./pd_sim.mat"
 filename_export_csv =  "./pd_sim.csv"
 sun = [1,1,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 = 1000
 numperiods = 8
+sensors = np.matrix([
+                    [ 1,0,0],
+                    [ -1,0,0],
+                    [0,1,0],
+                    [0,-1,0],
+                    [0,0,1],
+                    [0,0,-1]
+                    ])
 
 #calculate simduration so we see numperiods periods of the smalest angular speed
 simduration = numperiods*np.pi / min(filter(lambda x: x>0.0,[w_x,w_y,w_z])) if w_x or w_y or w_z else 1.0 #in seconds
 simtimestep = simduration / numsamples #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
+angle_earthdisc_visible = np.arcsin(float(earth_radius) / earth_distance )
+assert(angle_earthdisc_visible < np.pi/2.0)
+w_xyz = np.matrix([w_x,w_y,w_z])
+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)
+angle_sun_earth =  np.arccos(max(min(np.dot(sun,earth),1.0),-1.0))
+#print angle_sun_earth /np.pi*180, angle_earthdisc_visible /np.pi*180
 
 print "Simulation Input:"
 print "　 Duration:", simduration
-print "　 Step Size", simtimestep
-print "　 Angular Speed x-axis [rad/s]", w_x
-print "　 Angular Speed y-axis [rad/s]", w_y
-print "　 Angular Speed z-axis [rad/s]", w_z
+print "　 Step Size:", simtimestep
+print "　 Angular Speed x-axis [rad/s]:", w_x
+print "　 Angular Speed y-axis [rad/s]:", w_y
+print "　 Angular Speed z-axis [rad/s]:", w_z
+print "   Earth " + ("Light intensity: %f" % earth_reflection_intensity if earth_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 ""
 
 
-simtimeline = np.arange(0,simduration,simtimestep)
 
-w_xyz = np.matrix([w_x,w_y,w_z])
-
-angular_speed_sum_vectors = ( w_xyz.transpose() * simtimeline ).transpose() % (2 *np.pi)
-
-sensors = np.matrix([
-                    [ 1,0,0],
-                    [ -1,0,0],
-                    [0,1,0],
-                    [0,-1,0],
-                    [0,0,1],
-                    [0,0,-1]
-                    ])
-sun = np.matrix(sun / np.linalg.norm(sun))
 
 def getRotationMatrix(angles):
     x_angle = angles[0,0]
@@ -72,13 +86,36 @@ def rotateVector(vec_to_rotate, angles_vector):
 
 def sensor_model_rel_sensitivity_cos(angle):
     if angle < np.pi/2.0:
-        return np.cos(angle[0,0])
+        return np.cos(angle)
+    else:
+        return 0.0
+
+sensor_model = sensor_model_rel_sensitivity_cos
+
+def sun_earth_shadow(angle):
+    global angle_sun_earth, angle_earthdisc_visible, sensor_model
+    if angle_sun_earth < angle_earthdisc_visible:
+        return 0.0
     else:
+        return sensor_model(angle)
+
+sun_model = sun_earth_shadow
+
+def earthdisc_light_rel_intensity(angle):
+    global angle_sun_earth, angle_earthdisc_visible, sensor_model
+    if angle_sun_earth < angle_earthdisc_visible:
         return 0.0
+    #todo: factor in angle of actually lightened surface depending on angle_sun_earth and the angle between sat orbit and the sun
+    elif angle < angle_earthdisc_visible:
+        return 1.0
+    else:
+        return sensor_model(angle - angle_earthdisc_visible)
+
+earth_model = earthdisc_light_rel_intensity
 
 def _quote(str):
     return '"'+ re.sub(r'(?!\\)"','\\"',str)  +'"'
-        
+
 def exportCSV(dst_file_path, names, tline, smatrix):
     with open(dst_file_path, "w") as fh:
         fh.truncate()
@@ -86,24 +123,27 @@ def exportCSV(dst_file_path, names, tline, smatrix):
         for row in np.vstack((tline,smatrix.T)).T:
             fh.write(";".join(["%f" % x for x in row.T]) +"\n")
 
-phicostrace=[]
 senstrace=[]
+suntrace=[]
+earthtrace=[]
 for cur_rot in angular_speed_sum_vectors:
-    sens1=rotateVector(sensors[0],cur_rot)
     rot_sensors=rotateVector(sensors,cur_rot)
     #sinphi=np.cross(sens1,sun) # /1 / 1
-    cosphi=np.dot(rot_sensors,sun.transpose())
-    #print np.arccos(sincos)
-    #phi=np.sqrt(np.square(sinphi).sum())*float(np.sign(sinphi).prod())
-    #phisintrace.append(np.arcsin(np.sqrt(np.sum(np.square(sinphi)))))
-    angle=np.arccos(cosphi)
-    #print angle
-    #print [sensor_rel_sensitivity_cos(a) for a in angle]
-    phicostrace.append(angle)
-    senstrace.append([sensor_model_rel_sensitivity_cos(a) for a in angle])
+    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]
+    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])
+    senstrace.append(light_intensity_per_sensor)
+
+    suntrace.append(angles_to_sun)
+    earthtrace.append(angles_to_earth)
 
 senstrace = np.matrix(senstrace)
-legend_names=["Sensor%d %s" % (x+1, sensors[x]) for x in range(0,6)]
+suntrace = np.matrix(suntrace)
+earthtrace = np.matrix(earthtrace)
+legend_names=["Sensor%d %s" % (x+1, sensors[x].flatten()) for x in range(0,6)]
 
 print "Saving to " + filename_export_mat
 scipy.io.savemat(filename_export_mat, mdict={"legend":legend_names, "senstrace":senstrace,"simtimeline":simtimeline,"angular_speed_xyz":w_xyz}, oned_as="column")
@@ -141,10 +181,13 @@ spncols=2
 spnrows=3
 for s in range(0,6):
   pl.subplot(spnrows,spncols,s+1)
-  pl.plot(simtimeline, senstrace[:,s])
+  pl.plot(simtimeline, senstrace[:,s],"g")
+  pl.plot(simtimeline, 1/(2*np.pi)*suntrace[:,s],"y")
+  pl.plot(simtimeline, 1/(2*np.pi)*earthtrace[:,s],"b")
   pl.xticks(xlabelpos,xlabelvalues,rotation=0)
   pl.yticks(ylabelpos, ylabelvalues)
   pl.ylabel("[Light Intensity]")
   pl.title(legend_names[s])
+  pl.legend(["Sensor Value","Sun Angle", "Earth Angle"])
 
 pl.show()