mambo5.py

'''Author: Kristofer DrozdDate: November 20, 2015Description:This python script performs calculations to obtain positional parameters (azi & zen) of the spacecraft and sun at the exact moment a picture is taken of the landmark being analyzed.The following text files must be in the same directory in which the code is being run:LMRKNAMES.txt: name of the landmark being analyzedNUMBERPIC.txt: # of pictures taken of landmarkRESOLUTION.txt: resolution of each picturePICTIMES.txt: UTC time of each pictureLAT.txt: latitude of landmarkLON.txt: west longitude of landmarkSCOBJ1.txt: x component of space craft to object center vectors (BF frame)SCOBJ2.txt: y component of space craft to object center vectors (BF frame)SCOBJ3.txt: z component of space craft to object center vectors (BF frame)SZ1.txt: x component of object center to sun unit vectors (BF frame)SZ2.txt: y component of object center to sun unit vectors (BF frame)SZ3.txt: z component of object center to sun unit vectors (BF frame)VLM1.txt: x component of object center to landmark vectors (BF frame)VLM2.txt: y component of object center to landmark vectors (BF frame)VLM3.txt: z component of object center to landmark vectors (BF frame)Output:A text file of the name <landmark name>_viewing.txt is created in the directory the code is run'''import numpy as npimport matplotlib.pyplot as pltimport matplotlib.patches as mpatchesdef file_len(fname): ''' This function counts the number of lines in a text file Parameters fname: name of the text file Returns i + 1: number of lines in the text file ''' with open(fname) as f: for i, l in enumerate(f): pass return i + 1def SCPOS(SCOBJ1, SCOBJ2, SCOBJ3): ''' This function creates the object's center to spacecraft vector in body fixed frame. Since the SCOBJ compoents are from the spacecraft to the object's center, the components are multiplied by -1. Parameters SCOBJ1: x component of SCOBJ SCOBJ2: y component of SCOBJ SCOBGJ: z component of SCOBJ Returns SCOBJ: vector ''' SCPOS = np.array([-SCOBJ1, -SCOBJ2, -SCOBJ3]) return SCPOSdef SUNPOS(SZ1, SZ2, SZ3): ''' This function creates the object center to sun vector in body fixed frame. Since SZ is a unit vector it is multuplied by 1 AU. Parameters SZ1: x component of SZ SZ2: y component of SZ SZ3: z component of SZ Returns SCOBJ: vector ''' SZ_au = np.multiply(np.array([ SZ1, SZ2, SZ3]), 1.496e8) return SZ_audef BEN2LMRK(VLM1, VLM2, VLM3): ''' This function creates the object center to landmark vector in body fixed frame. Parameters VLM1: x component of BEN2LMRK VLM2: y component of BEN2LMRK VLM3: z component of BEN2LMRK Returns yomama: vector ''' yomama = np.array([VLM1, VLM2, VLM3], float) return yomama def BF2SEZ_tm(lon,lat): ''' This function creates the Body Fixed frame to SEZ fram transformation matrix. Parameters lon: The east longitude of the landamrk (degrees) lat: The latitude of the landamrk (degreesP) Returns Qxx_mat: The transformation matrix ''' lat_r = np.radians(lat) lon_r = np.radians(lon) Qxx_mat = np.array([[ np.sin(lat_r)*np.cos(lon_r), np.sin(lat_r)*np.sin(lon_r), -np.cos(lat_r)], [-np.sin(lon_r), np.cos(lon_r), 0], [np.cos(lat_r)*np.cos(lon_r), np.cos(lat_r)*np.sin(lon_r),np.sin(lat_r)]]) return Qxx_mat#==========================================================def main(): ''' The various input files are read in, so that their data can be used to calculate the azimuth and zenith angles of the spacecraft and sun with respect to a topcentric NEZ frame centered on the landmark during the exact moment a picture is taken of said landmark. A text file is then created listing the angles, picture resolution, and picture UTC time. ''' ''' Reading in TXT Files ''' fid1 = open('LAT.txt') fid2 = open('LON.txt') fid3 = open('NUMBERPIC.txt') fid4 = open('PICTIMES.txt') fid5 = open('lmrkname.txt') fid6 = open('RESOLUTION.txt') fid7 = open('SCOBJ1.txt') fid8 = open('SCOBJ2.txt') fid9 = open('SCOBJ3.txt') fid10 = open('SZ1.txt') fid11 = open('SZ2.txt') fid12 = open('SZ3.txt') fid13 = open('VLM1.txt') fid14 = open('VLM2.txt') fid15 = open('VLM3.txt') ''' Getting latitude and longitude of landmark ''' lat = float(fid1.readline().rstrip()) lon = 360 - float(fid2.readline().rstrip()) ''' Getting object center to landmark vector components and then vector ''' VLM1 = fid13.readline().rstrip() VLM1 = float(VLM1.replace('D', 'E')) VLM2 = fid14.readline().rstrip() VLM2 = float(VLM2.replace('D', 'E')) VLM3 = fid15.readline().rstrip() VLM3 = float(VLM3.replace('D', 'E')) lmrk_vec = BEN2LMRK(VLM1, VLM2, VLM3) ''' Getting number of pictures and landmark name ''' number_pics = int(fid3.readline().rstrip()) lmrk_name = fid5.readline().rstrip() ''' Formting outputs & creating text file ''' file = open(lmrk_name+"_viewing.txt", "w") file.write( "\n") file.write( lmrk_name) file.write( "\n") file.write( "lat = %f, Elon = %f\n" % (lat, lon)) file.write( "# of pictures taken = %f\n" % (number_pics)) file.write("\n") file.write( "| UTC | res | sun_zen | sun_azi | sc_zen | sc_azi |\n") file.write( "|--------------------------------------------------------------------------------------------|\n") ''' for loop ''' for j in range(0, number_pics): ''' Getting the SCOBJ vector components one line at a time ''' SCOBJ1 = fid7.readline().rstrip() SCOBJ1 = float(SCOBJ1.replace('D', 'E')) SCOBJ2 = fid8.readline().rstrip() SCOBJ2 = float(SCOBJ2.replace('D', 'E')) SCOBJ3 = fid9.readline().rstrip() SCOBJ3 = float(SCOBJ3.replace('D', 'E')) ''' Getting the SZ vector components one line at a time ''' SZ1 = fid10.readline().rstrip() SZ1 = float(SZ1.replace('D', 'E')) SZ2 = fid11.readline().rstrip() SZ2 = float(SZ2.replace('D', 'E')) SZ3 = fid12.readline().rstrip() SZ3 = float(SZ3.replace('D', 'E')) ''' Getting picture resolution one line at a time ''' res = float(fid6.readline().rstrip()) ''' Getting UTC time of each picture one line at a time ''' time = fid4.readline().rstrip() ''' SC Calculations lmrk2sc_vec: landmark to sc vector in body fixed frame n_sc: landmark to sc vector in topographical SEZ frame elvation_sc: elevation angle of n_sc zenith_sc: zenith angle of n_sc azimuth_sc: azimuth angle of n_sc, but manipulated so it is in NEZ frame opposed to SEZ ''' lmrk2sc_vec = np.subtract(SCPOS(SCOBJ1, SCOBJ2, SCOBJ3),lmrk_vec) n_sc = np.dot(BF2SEZ_tm(lon,lat),lmrk2sc_vec) elevation_sc = np.degrees(np.arcsin(np.true_divide(n_sc[2],np.linalg.norm(n_sc)))) zenith_sc = 90-elevation_sc azimuth_sc = np.degrees(np.arctan2(n_sc[1],-n_sc[0])) if azimuth_sc < 0: azimuth_sc = 360+azimuth_sc else: azimuth_sc = azimuth_sc ''' Sun Calculations lmrk2sun_vec: landmark to sun vector in body fixed frame n_sun: landmark to sun vector in topographical SEZ frame elvation_sun: elevation angle of n_sun zenith_sun: zenith angle of n_sun azimuth_sun: azimuth angle of n_sun, but manipulated so it is in NEZ frame opposed to SEZ ''' lmrk2sun_vec = np.subtract(SUNPOS(SZ1, SZ2, SZ3),lmrk_vec) n_sun = np.dot(BF2SEZ_tm(lon,lat),lmrk2sun_vec) elevation_sun = np.degrees(np.arcsin(np.true_divide(n_sun[2],np.linalg.norm(n_sun)))) zenith_sun = 90-elevation_sun azimuth_sun = np.degrees(np.arctan2(n_sun[1],-n_sun[0])) if azimuth_sun < 0: azimuth_sun = 360+azimuth_sun else: azimuth_sun = azimuth_sun ''' Continued formatting of text file ''' file.write("| ") file.write(time) file.write(" ") file.write( "| %07.3f | %+08.3f | %08.3f | %+08.3f | %08.3f \n" % (res, zenith_sun, azimuth_sun, zenith_sc, azimuth_sc)) file.close() with open(lmrk_name+"_viewing.txt") as ifh: arr = np.loadtxt(ifh, usecols = (1,2,3,4,5), dtype = float, delimiter = " | ", skiprows = 7) theta = np.linspace(-np.pi, np.pi, 100) azi_sun_1 = arr[:,2] zen_sun_1 = arr[:,1] azi_sc_1 = arr[:,4] zen_sc_1 = arr[:,3] res = arr[:,0] ''' Sun plot ''' plt.figure(1) ax1 = plt.subplot(111, polar = True) ax1.set_theta_zero_location("N") ax1.set_theta_direction(-1) plt.grid(True) ax1.set_rgrids([15,30,45,60,75,90], angle = 60, fontsize = 10) ax1.set_rlim(0, 90) ax1.set_thetagrids([0, 45, 90, 135, 180, 225, 270, 315], frac = 1.08, fontsize = 10) for k in range(0, number_pics): if res[k] >=1: ax1.plot(np.radians(azi_sun_1[k]), zen_sun_1[k], 'or') elif 1 > res[k] and res[k] >= .5: ax1.plot(np.radians(azi_sun_1[k]), zen_sun_1[k], 'og') elif .5 > res[k] and res[k] >= .25: ax1.plot(np.radians(azi_sun_1[k]), zen_sun_1[k], 'ob') elif .25 > res[k] and res[k] >= .1: ax1.plot(np.radians(azi_sun_1[k]), zen_sun_1[k], 'oy') elif .1 > res[k] and res[k] >= .05: ax1.plot(np.radians(azi_sun_1[k]), zen_sun_1[k], 'om') else: ax1.plot(np.radians(azi_sun_1[k]), zen_sun_1[k], 'oc') mot1 = mpatches.Patch( color = 'r', label = 'Res $\geq$ 50 cm') mot2 = mpatches.Patch( color = 'g', label = '100 cm > Res $\geq$ 50 cm') mot3 = mpatches.Patch( color = 'b', label = '50 cm > Res $\geq$ 25 cm') mot4 = mpatches.Patch( color = 'y', label = '25 cm > Res $\geq$ 10 cm') mot5 = mpatches.Patch( color = 'm', label = '10 cm > Res $\geq$ 5 cm') mot6 = mpatches.Patch( color = 'c', label = 'Res < 5 cm') plt.legend( handles = [mot1, mot2, mot3, mot4, mot5, mot6], fontsize = 8, loc = 2, bbox_to_anchor = (.93, 1.1)) figure_title1 = lmrk_name+' Sun Azimuth vs. Zenith' plt.text(.5, 1.08,figure_title1, horizontalalignment = 'center', fontsize = 15, transform = ax1.transAxes) plt.savefig( lmrk_name+'_sun_plot.png') plt.clf() ''' SC plot ''' plt.figure(2) ax2 = plt.subplot(111, polar = True) ax2.set_theta_zero_location("N") ax2.set_theta_direction(-1) plt.grid(True) ax2.set_rgrids([15,30,45,60,75,90], angle = 60, fontsize = 10) ax2.set_rlim(0, 90) ax2.set_thetagrids([0, 45, 90, 135, 180, 225, 270, 315], frac = 1.08, fontsize = 10) for k in range(0, number_pics): if res[k] >=1: ax2.plot(np.radians(azi_sc_1[k]), zen_sc_1[k], 'or') elif 1 > res[k] and res[k] >= .5: ax2.plot(np.radians(azi_sc_1[k]), zen_sc_1[k], 'og') elif .5 > res[k] and res[k] >= .25: ax2.plot(np.radians(azi_sc_1[k]), zen_sc_1[k], 'ob') elif .25 > res[k] and res[k] >= .1: ax2.plot(np.radians(azi_sc_1[k]), zen_sc_1[k], 'oy') elif .1 > res[k] and res[k] >= .05: ax2.plot(np.radians(azi_sc_1[k]), zen_sc_1[k], 'om') else: ax2.plot(np.radians(azi_sc_1[k]), zen_sc_1[k], 'oc') dot1 = mpatches.Patch( color = 'r', label = 'Res $\geq$ 50 cm') dot2 = mpatches.Patch( color = 'g', label = '100 cm > Res $\geq$ 50 cm') dot3 = mpatches.Patch( color = 'b', label = '50 cm > Res $\geq$ 25 cm') dot4 = mpatches.Patch( color = 'y', label = '25 cm > Res $\geq$ 10 cm') dot5 = mpatches.Patch( color = 'm', label = '10 cm > Res $\geq$ 5 cm') dot6 = mpatches.Patch( color = 'c', label = 'Res < 5 cm') plt.legend( handles = [dot1, dot2, dot3, dot4, dot5, dot6], fontsize = 8, loc = 2, bbox_to_anchor = (.93, 1.1)) figure_title2 = lmrk_name+' SC Azimuth vs. Zenith' plt.text(.5, 1.08,figure_title2, horizontalalignment = 'center', fontsize = 15, transform = ax2.transAxes) plt.savefig(lmrk_name+'_sc_plot.png') plt.clf() #print(n_sc) #print(lmrk2sc_vec) #print(elevation_sc) #print(elevation_sun) #print(BF2SEZ_tm(lon,lat))if __name__ == '__main__': main()