stud

study spacejunk
Log | Files | Refs | LICENSE

task2_1.py (2823B) - Raw

      1 #!/usr/bin/python3
      2 
      3 import sys
      4 import csv
      5 
      6 from math import radians, degrees, tan, atan, pi, log, sin, cos
      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 rphi_1, rphi_2 = radians(phi_1), radians(phi_2)
     40 rlambda_1, rlambda_2 = radians(lambda_1), radians(lambda_2)
     41 phil = round((phi_p+phi_s)/2)
     42 nphi = int((phi_s-phi_p)/dphi)+1
     43 nlambda = int((lambda_r-lambda_v)/dlambda)+1
     44 
     45 # label orientations
     46 W, E, N, S = (-25, -5), (10, -5), (-5, 10), (-5, -20)
     47 SW, NE = (-10, -10), (5, 5)
     48 
     49 kr = {}
     50 with open("krasovskio.csv") as f:
     51     for row in csv.DictReader(f):
     52         kr[float(row['phi'])] = row
     53 betamm = float(kr[phil]["r"]) * 1000 / M
     54 
     55 
     56 def yx(lat, lon):
     57     # lat - phi in degrees
     58     # lon - lambda in degrees
     59     phi = radians(lat)
     60     U = tan(pi/4 + phi/2)
     61     xmm = betamm * log(U)
     62     ymm = betamm * radians(lon)
     63     return (ymm, xmm)
     64 
     65 
     66 points = []
     67 for i in range(nphi):
     68     phid = phi_p + i*dphi
     69     on_y = []
     70     for j in range(nlambda):
     71         lambdad = lambda_v + j*dlambda
     72         on_y.append(yx(phid, lambdad))
     73     points.append(sorted(on_y))
     74 
     75 
     76 fig, ax = plt.subplots()
     77 ax.set_aspect('equal')
     78 ax.axis("off")
     79 
     80 
     81 # abscises
     82 for i in range(nphi):
     83     row = [points[i][j] for j in range(nlambda)]
     84     ax.plot(*(LineString(row).xy), color="xkcd:black", linewidth=.5)
     85     annotate(ax, c(phi_p+i*dphi), row[0], W)
     86     annotate(ax, c(phi_p+i*dphi), row[-1], E)
     87 
     88 # ordinates
     89 for i in range(nlambda):
     90     col = [points[j][i] for j in range(nphi)]
     91     ax.plot(*(LineString(col).xy), color="xkcd:black", linewidth=.5)
     92     annotate(ax, c(lambda_v+i*dlambda), col[0], S)
     93     annotate(ax, c(lambda_v+i*dlambda), col[-1], N)
     94 
     95 # loksodroma
     96 A = yx(phi_1, lambda_1)
     97 B = yx(phi_2, lambda_2)
     98 loksodroma = ((A, B))
     99 ax.plot(*(LineString(loksodroma).xy), color="xkcd:black", linewidth=.5)
    100 
    101 # ortodroma
    102 ctgu = ctg(rphi_1)*tan(rphi_2)*cosec(rlambda_2-rlambda_1)-ctg(rlambda_2-rlambda_1)
    103 u = arccot(ctgu)
    104 ortodroma = []
    105 for lambdad in lambda_range:
    106     phi_ort = atan(tan(rphi_1)*cosec(u)*sin(u-rlambda_1+radians(lambdad)))
    107     ortodroma.append(yx(round(degrees(phi_ort)*2)/2, lambdad))
    108 ax.plot(*(LineString(ortodroma).xy), color="xkcd:black", linewidth=.5)
    109 
    110 annotate(ax, "A", A, SW)
    111 annotate(ax, "B", B, NE)
    112 
    113 if __name__ == '__main__':
    114     if len(sys.argv) == 2:
    115         plt.savefig(sys.argv[1], bbox_inches='tight')
    116         print("Saved %s" % sys.argv[1])
    117     else:
    118         plt.show()