*-----------------------------------------------------------------------
*     MAP PROJECTION (PTOLEMAIC 2ND)                2007-11-19 E. TOYODA
*-----------------------------------------------------------------------
*     Copyright (C) 2000-2016 GFD Dennou Club. All rights reserved.
*-----------------------------------------------------------------------
      SUBROUTINE MPFPT2(XLON, YLAT, X, Y)

      PARAMETER (EPSL = 1.0E-5)
      PARAMETER (PI = ATAN(1.0) * 4.0)
      PARAMETER (HALFPI = ATAN(1.0) * 2.0)
      PARAMETER (THULE = 63.*PI/180.0)
      PARAMETER (SYENE = (23.+50./60.)*PI/180.0)
      REAL MEROE
      PARAMETER (MEROE = -(16.+25./60.)*PI/180.0)

      EXTERNAL  XMPLON

        XX = XMPLON(XLON)
        R1 = THULE + HALFPI
        TH1 = COS(THULE) * XX / (PI - THULE)
        X1 = R1 * SIN(TH1)
        Y1 = PI - R1 * COS(TH1)
        print '(a,"(",2f12.7,")")', 'x1,y1=', x1, y1
        R2 = SYENE + HALFPI
        TH2 = COS(SYENE) * XX / (PI - SYENE)
        X2 = R2 * SIN(TH2)
        Y2 = PI - R2 * COS(TH2)
        print '(a,"(",2f12.7,")")', 'x2,y2=', x2, y2
        R3 = MEROE + HALFPI
        TH3 = COS(MEROE) * XX / (PI - MEROE)
        X3 = R3 * SIN(TH3)
        Y3 = PI - R3 * COS(TH3)
        print '(a,"(",2f12.7,")")', 'x3,y3=', x3, y3
        DET = (Y1 - Y2) * (X2 - X3) - (X2 - X1) * (Y3 - Y2)
        print '(a,f12.7)', 'det:', det
        C1 = 0.5 * ((X2 - X3)*(X3 - X1) + (Y2 - Y3)*(Y3 - Y1)) / DET
        C2 = 0.5 * ((X1 - X2)*(X3 - X1) + (Y1 - Y2)*(Y3 - Y1)) / DET
        print '(a,"(",2f12.7,")")', 'c1,c2=', c1, c2
        XC1 = 0.5 * (X1 + X2) + C1 * (Y1 - Y2)
        YC1 = 0.5 * (Y1 + Y2) + C1 * (X2 - X1)
        XC2 = 0.5 * (X2 + X3) + C2 * (Y2 - Y3)
        YC2 = 0.5 * (Y2 + Y3) + C2 * (X3 - X2)
        XC = 0.5 * (XC1 + XC2)
        YC = 0.5 * (YC1 + YC2)
        print '(a,"(",2f12.7,")")', 'center=', xc, yc
        R = SQRT((X2 - XC) ** 2 + (Y2 - YC) ** 2)
        R0 = PI - YLAT
        RHO = SQRT(XC ** 2 + (PI - YC) ** 2)
        print '(a,3f12.7)', 'r,r0,rho=', r, r0, RHO
        IF (RHO .GT. R0 + R) THEN
          CALL GLRGET('RUNDEF',RNA)
          X = RNA
          Y = RNA
          RETURN
        ENDIF
        ALPHA = ACOS((-R**2 + RHO**2 + R0**2) / (2.0 * RHO * R0))
        IF (ATAN2(X2, -Y2) .LT. ATAN2(XC, -YC)) THEN
          ALPHA = -ALPHA
        ENDIF
        print '(a,f12.7)', 'alpha=', alpha
        R0RHO = R0 / RHO
        X = R0RHO * (COS(ALPHA) * XC - SIN(ALPHA) * (YC - PI))
        Y = R0RHO * (SIN(ALPHA) * XC + COS(ALPHA) * (YC - PI)) + PI
      RETURN
*-----------------------------------------------------------------------
      ENTRY MPIPT2(X, Y, XLON, YLAT)

        XLON = X
        YLAT = Y
        IF (ABS(XLON) .LE. PI) RETURN

        CALL GLRGET('RUNDEF',RNA)
        XLON = RNA
        YLAT = RNA
      RETURN
      END
