task2_2.py (3608B) - Raw
1 #!/usr/bin/python3 2 3 import sys 4 import csv 5 from math import radians, degrees, tan, atan, sin, cos 6 7 from shapely.geometry import LineString 8 import matplotlib.pyplot as plt 9 10 from consts import ( 11 phi_p, phi_s, dphi, 12 phi_1, phi_2, 13 lambda_1, lambda_2, lambda_range, 14 lambda_v, lambda_r, dlambda, 15 M 16 ) 17 18 19 def c(x): 20 return "{}°".format(x) 21 22 23 def annotate(ax, text, point, heading): 24 ax.annotate(text, point, textcoords="offset points", xytext=heading) 25 26 27 def ctg(x): 28 return cos(x)/sin(x) 29 30 31 def arccot(x): 32 return atan(1/x) 33 34 35 def cosec(x): 36 return 1/sin(x) 37 38 39 phi_loks = 27.308 # loksodromos platuma 12 ilgumoje interpoliavus 40 41 rphi_1, rphi_2 = radians(phi_1), radians(phi_2) 42 rlambda_1, rlambda_2 = radians(lambda_1), radians(lambda_2) 43 phil = round((phi_p+phi_s)/2) 44 nphi = int((phi_s-phi_p)/dphi)+1 45 nlambda = int((lambda_r-lambda_v)/dlambda)+1 46 midlambda = int(lambda_r+lambda_v)/2 47 midnlambda = int((lambda_r-midlambda)/dlambda)+1 48 49 # label orientations 50 W, E, N, S = (-25, -5), (10, -5), (-5, 10), (-5, -20) 51 SW, NE = (-10, -10), (5, 5) 52 53 kr = {} 54 with open("krasovskio.csv") as f: 55 for row in csv.DictReader(f): 56 kr[float(row['phi'])] = {k: float(v) for k, v in row.items()} 57 58 alpha = (kr[phi_1]["lgr"]-kr[phi_2]["lgr"])/(kr[phi_2]["lgU"]-kr[phi_1]["lgU"]) 59 Ualpha = (10**kr[phi_p]["lgU"])**alpha 60 U1alpha = (10**kr[phi_1]["lgU"])**alpha 61 U2alpha = (10**kr[phi_2]["lgU"])**alpha 62 C1 = (kr[phi_1]["r"]*U1alpha)/alpha 63 C2 = (kr[phi_2]["r"]*U2alpha)/alpha 64 if abs(C1 - C2) / C1 > 1e-6: 65 raise ValueError("too large error between C1 and C2") 66 Cmm = C1 * 1000 / M 67 qmm = Cmm/Ualpha 68 69 70 def yx(lat, lon): 71 # lat - phi in degrees 72 # lon - lambda in degrees 73 lgU = kr[round(lat*2)/2.]["lgU"] 74 Ualpha = (10**lgU)**alpha 75 pmm = Cmm/Ualpha 76 delta = alpha * lon 77 xmm = qmm-pmm*cos(radians(delta)) 78 ymm = pmm*sin(radians(delta)) 79 return (ymm, xmm) 80 81 82 points = [] 83 for i in range(nphi): 84 phid = phi_p + i*dphi 85 on_y = [] 86 for j in range(midnlambda): 87 lambdad = j*dlambda 88 ymm, xmm = yx(phid, lambdad) 89 on_y.append((ymm, xmm)) 90 if j > 0: 91 on_y.append((-ymm, xmm)) 92 points.append(sorted(on_y)) 93 94 95 fig, ax = plt.subplots() 96 ax.set_aspect('equal') 97 ax.axis("off") 98 99 100 # abscises 101 for i in range(nphi): 102 row = [points[i][j] for j in range(nlambda)] 103 ax.plot(*(LineString(row).xy), color="xkcd:black", linewidth=.5) 104 annotate(ax, c(phi_p+i*dphi), row[0], W) 105 annotate(ax, c(phi_p+i*dphi), row[-1], E) 106 107 # ordinates 108 for i in range(nlambda): 109 col = [points[j][i] for j in range(nphi)] 110 ax.plot(*(LineString(col).xy), color="xkcd:black", linewidth=.5) 111 annotate(ax, c(lambda_v+i*dlambda), col[0], S) 112 annotate(ax, c(lambda_v+i*dlambda), col[-1], N) 113 114 # loksodroma 115 rmidlambda = radians(midlambda) 116 A = yx(phi_1, lambda_1-midlambda) 117 MidLoks = yx(phi_loks, 0) 118 B = yx(phi_2, lambda_2-midlambda) 119 loksodroma = ((A, MidLoks, B)) 120 ax.plot(*(LineString(loksodroma).xy), color="xkcd:black", linewidth=.5) 121 122 # ortodroma 123 ctgu = ctg(rphi_1)*tan(rphi_2)*cosec(rlambda_2-rlambda_1)-ctg(rlambda_2-rlambda_1) 124 u = arccot(ctgu) 125 ortodroma = [] 126 for lambdad in lambda_range: 127 phi_ort = atan(tan(rphi_1)*cosec(u)*sin(u-rlambda_1+radians(lambdad))) 128 ortodroma.append(yx(round(degrees(phi_ort)*2)/2, lambdad-midlambda)) 129 ax.plot(*(LineString(ortodroma).xy), color="xkcd:black", linewidth=.5) 130 131 annotate(ax, "A", A, SW) 132 annotate(ax, "B", B, NE) 133 134 if __name__ == '__main__': 135 if len(sys.argv) == 2: 136 plt.savefig(sys.argv[1], bbox_inches='tight') 137 print("Saved %s" % sys.argv[1]) 138 else: 139 plt.show()