1*-----------------------------------------------------------------------
2*     MAP PROJECTION (PTOLEMAIC 2ND)                2007-11-19 E. TOYODA
3*-----------------------------------------------------------------------
4*     Copyright (C) 2000-2007 GFD Dennou Club. All rights reserved.
5*-----------------------------------------------------------------------
6      SUBROUTINE MPFPT2(XLON, YLAT, X, Y)
7
8      PARAMETER (EPSL = 1.0E-5)
9      PARAMETER (PI = ATAN(1.0) * 4.0)
10      PARAMETER (HALFPI = ATAN(1.0) * 2.0)
11      PARAMETER (THULE = 63.*PI/180.0)
12      PARAMETER (SYENE = (23.+50./60.)*PI/180.0)
13      REAL MEROE
14      PARAMETER (MEROE = -(16.+25./60.)*PI/180.0)
15
16      EXTERNAL  XMPLON
17
18        XX = XMPLON(XLON)
19        R1 = THULE + HALFPI
20        TH1 = COS(THULE) * XX / (PI - THULE)
21        X1 = R1 * SIN(TH1)
22        Y1 = PI - R1 * COS(TH1)
23        print '(a,"(",2f12.7,")")', 'x1,y1=', x1, y1
24        R2 = SYENE + HALFPI
25        TH2 = COS(SYENE) * XX / (PI - SYENE)
26        X2 = R2 * SIN(TH2)
27        Y2 = PI - R2 * COS(TH2)
28        print '(a,"(",2f12.7,")")', 'x2,y2=', x2, y2
29        R3 = MEROE + HALFPI
30        TH3 = COS(MEROE) * XX / (PI - MEROE)
31        X3 = R3 * SIN(TH3)
32        Y3 = PI - R3 * COS(TH3)
33        print '(a,"(",2f12.7,")")', 'x3,y3=', x3, y3
34        DET = (Y1 - Y2) * (X2 - X3) - (X2 - X1) * (Y3 - Y2)
35        print '(a,f12.7)', 'det:', det
36        C1 = 0.5 * ((X2 - X3)*(X3 - X1) + (Y2 - Y3)*(Y3 - Y1)) / DET
37        C2 = 0.5 * ((X1 - X2)*(X3 - X1) + (Y1 - Y2)*(Y3 - Y1)) / DET
38        print '(a,"(",2f12.7,")")', 'c1,c2=', c1, c2
39        XC1 = 0.5 * (X1 + X2) + C1 * (Y1 - Y2)
40        YC1 = 0.5 * (Y1 + Y2) + C1 * (X2 - X1)
41        XC2 = 0.5 * (X2 + X3) + C2 * (Y2 - Y3)
42        YC2 = 0.5 * (Y2 + Y3) + C2 * (X3 - X2)
43        XC = 0.5 * (XC1 + XC2)
44        YC = 0.5 * (YC1 + YC2)
45        print '(a,"(",2f12.7,")")', 'center=', xc, yc
46        R = SQRT((X2 - XC) ** 2 + (Y2 - YC) ** 2)
47        R0 = PI - YLAT
48        RHO = SQRT(XC ** 2 + (PI - YC) ** 2)
49        print '(a,3f12.7)', 'r,r0,rho=', r, r0, RHO
50        IF (RHO .GT. R0 + R) THEN
51          CALL GLRGET('RUNDEF',RNA)
52          X = RNA
53          Y = RNA
54          RETURN
55        ENDIF
56        ALPHA = ACOS((-R**2 + RHO**2 + R0**2) / (2.0 * RHO * R0))
57        IF (ATAN2(X2, -Y2) .LT. ATAN2(XC, -YC)) THEN
58          ALPHA = -ALPHA
59        ENDIF
60        print '(a,f12.7)', 'alpha=', alpha
61        R0RHO = R0 / RHO
62        X = R0RHO * (COS(ALPHA) * XC - SIN(ALPHA) * (YC - PI))
63        Y = R0RHO * (SIN(ALPHA) * XC + COS(ALPHA) * (YC - PI)) + PI
64      RETURN
65*-----------------------------------------------------------------------
66      ENTRY MPIPT2(X, Y, XLON, YLAT)
67
68        XLON = X
69        YLAT = Y
70        IF (ABS(XLON) .LE. PI) RETURN
71
72        CALL GLRGET('RUNDEF',RNA)
73        XLON = RNA
74        YLAT = RNA
75      RETURN
76      END
77