From 4b3c473068f09e70fc37be33c76ba969dac43e6e Mon Sep 17 00:00:00 2001 From: Bernhard Tittelbach Date: Wed, 21 Mar 2012 16:02:42 +0000 Subject: added earthdisc to photodiode simulator git-svn-id: https://svn.spreadspace.org/mur.sat@324 7de4ea59-55d0-425e-a1af-a3118ea81d4c --- tools/photodiodesim/pd_sim.py | 111 +++++++++++++++++++++++++++++------------- 1 file changed, 77 insertions(+), 34 deletions(-) (limited to 'tools/photodiodesim/pd_sim.py') 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() -- cgit v1.2.3