GPLK map projection

import math

pi=3.14159265358979
degrees=pi/180

############################ ENTER DATA ############################
mapin=open('J:\EX-IM\W2H\Map projection\Atlas_of_of_Venus_(Venusian_Haven).V1.svg', 'rb')
mapout=open('J:\EX-IM\W2H\Map projection\Atlas_of_of_Venus_(Venusian_Haven).EasternHemisphere.svg', 'wb')
C5=-6*degrees # Point A longitude
C6=0*degrees # Point A latitude
D5=173*degrees # Point B longitude
D6=-0*degrees # Point B latitude
C13=0 # 1 for Lambert azimuthal, 0 for Gall-Peters cilindric
####################################################################

l4c="1234" # the last four characters
readcoordinates=False
watchit=False
numberstatus=0  # 0: not within a coordinate pair
                # 1: read the first number
                # 2: between the two coordinates
                # 3: read the second number
num1=''
num2=''
numchars=['0','1','2','3','4','5','6','7','8','9','.','-','+','e']

def roundit(sn=3.14159265358979):
    rnum=round(sn,5) # the number stands for the number of digits after the point
    position=`sn`.find('.')
    snum=`sn`[:position+6]
    return snum

C7=math.cos(C5)*math.cos(C6)
C8=math.sin(C5)*math.cos(C6)
C9=math.sin(C6)
D7=math.cos(D5)*math.cos(D6)
D8=math.sin(D5)*math.cos(D6)
D9=math.sin(D6)
E7=(C7+D7)/2
E8=(C8+D8)/2
E9=(C9+D9)/2
E10=math.sqrt(E7*E7+E8*E8+E9*E9)
F7=E7/E10
F8=E8/E10
F9=E9/E10
F5=math.atan2(F8,F7)
F6=math.asin(F9)

G3=math.cos(F5)*math.cos(F6)
G4=math.sin(F5)*math.cos(F6)
G5=math.sin(F6)
H3=-math.sin(F5)
H4=math.cos(F5)
H5=0
I3=-math.sin(F6)*math.cos(F5)
I4=-math.sin(F6)*math.sin(F5)
I5=math.cos(F6)

D16=C5
E16=C6
F16=math.cos(D16)*math.cos(E16)
G16=math.sin(D16)*math.cos(E16)
H16=math.sin(E16)
I16=F16*H4*I5+G16*H5*I3+H16*H3*I4-H16*H4*I3-H5*I4*F16-I5*G16*H3
J16=G3*G16*I5+G4*H16*I3+G5*F16*I4-G5*G16*I3-H16*I4*G3-I5*G4*F16
K16=G3*H4*H16+G4*H5*F16+G5*H3*G16-G5*H4*F16-H5*G16*G3-H16*G4*H3
L16=math.atan2(J16,I16)
M16=math.asin(K16)
P16=-math.atan2(J16,K16)*1-pi/2

G7=math.cos(F5)*math.cos(F6)
G8=math.sin(F5)*math.cos(F6)
G9=math.sin(F6)
H7=math.cos(P16)*H3+math.sin(P16)*I3
H8=math.cos(P16)*H4+math.sin(P16)*I4
H9=math.cos(P16)*H5+math.sin(P16)*I5
I7=math.cos(P16)*I3-math.sin(P16)*H3
I8=math.cos(P16)*I4-math.sin(P16)*H4
I9=math.cos(P16)*I5-math.sin(P16)*H5

def calculatecoordinates(D20,E20):
    F20=math.cos(D20)*math.cos(E20)
    G20=math.sin(D20)*math.cos(E20)
    H20=math.sin(E20)
    I20=F20*H8*I9+G20*H9*I7+H20*H7*I8-H20*H8*I7-H9*I8*F20-I9*G20*H7
    J20=G7*G20*I9+G8*H20*I7+G9*F20*I8-G9*G20*I7-H20*I8*G7-I9*G8*F20
    K20=G7*H8*H20+G8*H9*F20+G9*H7*G20-G9*H8*F20-H9*G20*G7-H20*G8*H7
    L20=math.atan2(J20,I20)
    M20=math.asin(K20)
    N20=L20
    O20=math.sin(M20)
    P20=math.atan2(J20,K20)
    Q20=math.sqrt(J20*J20+K20*K20+(1-I20)*(1-I20))
    R20=math.sin(P20)*Q20
    S20=math.cos(P20)*Q20
    T20=R20*C13+N20*(1-C13)
    U20=S20*C13+O20*(1-C13)
    W20=T20/degrees
    X20=U20/degrees
    Z20=math.atan2(X20,W20)+P16
    AA20=math.sqrt(W20*W20+X20*X20)
    AC20=math.cos(Z20)*AA20
    AD20=-math.sin(Z20)*AA20
    if watchit:
        print 'num1', num1
        print 'num2', num2
        print D20
        print E20
        print F20
        print G20
        print H20
        print I20
        print J20
        print K20
        print L20
        print M20
        print N20
        print O20
        print P20
        print Q20
        print R20
        print S20
        print T20
        print U20
        print W20
        print X20
        print Z20
        print AA20
        print AC20
        print AD20
    mapout.write(roundit(AC20)+","+roundit(AD20))

while True:
    try:
        ch=mapin.read(1)
        if ch=='': break
        if ch=='$': watchit=True
        l4c=l4c[1:]+ch
        if readcoordinates and (ch=='"'):
            readcoordinates=False
            if (l4c[2] in numchars) and (numberstatus==3):
                numberstatus=0
                calculatecoordinates((eval(num1)-180.0)*pi/180,(90.0-eval(num2))*pi/180)
                num1=""
                num2=""
        if l4c==' d="': readcoordinates=True
        if readcoordinates:
            if (numberstatus==0) and (ch in numchars): numberstatus=1
            if (numberstatus==1) and (not (ch in numchars)): numberstatus=2
            if (numberstatus==2) and (ch in numchars): numberstatus=3
            if (numberstatus==3) and (not (ch in numchars)):
                numberstatus=0
                calculatecoordinates((eval(num1)-180.0)*pi/180,(90.0-eval(num2))*pi/180)
                num1=""
                num2=""
            if numberstatus==1: num1=num1+ch
            if numberstatus==3: num2=num2+ch
            if numberstatus==0: mapout.write(ch)
        else:
            mapout.write(ch)
    except KeyboardInterrupt: break

mapout.close()
mapin.close()