C @(#)adxy.for	17.1.1.1 (ESO-DMD) 01/25/02 17:13:33
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
c----------------------------------------------------------------------
      SUBROUTINE ADXY(A0,D0,XEST,YEST,NITER,INEX,COSD0,SIND0,ALFA0,
     *NX,NY,KK,IXY,IDX,IDY,BX,BY,CX,CY,SCHMIDT)
c----------------------------------------------------------------------
      IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N)
C 
      PARAMETER   (MAXOBJ=100001)
C 
      DIMENSION A0(3),D0(3),D0A(3),AR(MAXOBJ,4),BX(9),BY(9),
     * XEST(2,2),YEST(2,2),X(10),KK(9),IDX(9),IDY(9)
      CHARACTER*1 SCHMIDT
C 
      PI = 4.0D0*DATAN(1.0D0)
      DOPI18 = PI/1.8D2
      DO0 = 0.0D0
      DO1 = 1.0D0
      DO60 = 6.0D1
      DO3600 = 3.6D3
      NITER = 0 
      INEX = 0
      ALFAS = (A0(1)+A0(2)/DO60+A0(3)/DO3600)*1.5D1*DOPI18
      RH = DO1

      DO 20 I = 1,3 
         IF (D0(I).LT.DO0) RH = -DO1 
   20    D0A(I) = DABS(D0(I))

      DELTS = RH*(D0A(1)+D0A(2)/DO60+D0A(3)/DO3600)*DOPI18
      IF (ALFA0.LT.PI/2D0.AND.ALFAS.GT.1.5D0*PI) ALFAS = ALFAS-2D0*PI 
      IF (ALFA0.GT.1.5D0*PI.AND.ALFAS.LT.PI/2D0) ALFAS = ALFAS+2D0*PI 
      COSRH = COSD0*DCOS(DELTS)*DCOS(ALFAS-ALFA0) + SIND0*DSIN(DELTS) 
      SINRH = DSQRT(DO1-COSRH**2) 
      SINPH = DCOS(DELTS)*DSIN(ALFAS-ALFA0)/SINRH 
      COSPH = (DSIN(DELTS)-SIND0*COSRH)/(COSD0*SINRH)
      IF (SCHMIDT.EQ.'Y') THEN 
          RH = DATAN(SINRH/COSRH)
      ELSE
          RH = SINRH/COSRH
      ENDIF 

      AKSI = RH*SINPH/DOPI18
      ANU = RH*COSPH/DOPI18 
      XEST(1,1) = ((AKSI-CX)*BY(2)-(ANU-CY)*BX(2))/ 
     1(BX(1)*BY(2)-BX(2)*BY(1)) 
      YEST(1,1) = ((AKSI-CX)*BY(1)-(ANU-CY)*BX(1))/ 
     1(BX(2)*BY(1)-BX(1)*BY(2)) 


   30 XEST(2,1) = XEST(1,1)*1.0D3 
      YEST(2,1) = YEST(1,1)*1.0D3 

      IF (DABS(XEST(2,1)).GT.2.0D5.OR.DABS(YEST(2,1)).GT.2.0D5) 
     >   then
         WRITE(6,9050) XEST(2,1),YEST(2,1) 
 9050    FORMAT(/,1X,'COMPUTED (X,Y):',2(F12.2),/, 
     >          1X,'POSITION OUTSIDE PLATE AREA. SORRY')
         INEX = 1
         return
         endif

   60 NITER = NITER + 1 
   70 AR(1,1) = XEST(1,1) 
      AR(1,2) = YEST(1,1)
      AR(1,3) = DO0 
      AR(1,4) = DO0 
      I16 = 1 
      N = NX+1
      IXY = 1 

      CALL RAR(I16,0,N,KK,IXY,IDX,IDY,X,AR) 

      L = NX
      V = DO0 

      DO 80 I17 = 1,L 
   80   V = V + BX(I17)*X(I17)
      AKSI1 = V + CX
      N = NY + 1
      IXY = 2 

      CALL RAR(I16,0,N,KK,IXY,IDX,IDY,X,AR) 

      L = NY
      V = DO0 

      DO 90 I17 = 1,L 
   90   V = V + BY(I17)*X(I17)

      ANU1 = V + CY 
      DKSI = AKSI1 - AKSI 
      DNU = ANU1 - ANU
      DDIST = DSQRT(DKSI**2 + DNU**2) 
      DCRIT = 1.0D-6
      IF (DDIST.LE.DCRIT ) return

      N = NX + 1
      IXY = 1 
      CALL RAR(1,1,N,KK,IXY,IDX,IDY,X,AR) 
      L = NX
      V = DO0 

      DO 110 I17 = 1,L
  110   V = V + BX(I17)*X(I17)

      DKDX = V
      N = NX + 1
      IXY = 1 
      CALL RAR(1,2,N,KK,IXY,IDX,IDY,X,AR) 
      L = NX
      V = DO0 

      DO 120 I17 = 1,L
  120   V = V + BX(I17)*X(I17)

      DKDY = V
      N = NY + 1
      IXY = 2 
      CALL RAR(1,1,N,KK,IXY,IDX,IDY,X,AR) 
      L = NY
      V = DO0 

      DO 130 I17 = 1,L
  130 V = V + BY(I17)*X(I17)

      DNDX = V
      N = NY + 1
      IXY = 2 
      CALL RAR(1,2,N,KK,IXY,IDX,IDY,X,AR) 
      L = NY
      V = DO0 

      DO 140 I17 = 1,L
  140   V = V + BY(I17)*X(I17)

      DNDY = V
      YEST(1,2) = (DKSI*DNDX-DNU*DKDX)/(DNDX*DKDY-DKDX*DNDY)
      XEST(1,2) = (DKSI*DNDY-DNU*DKDY)/(DKDX*DNDY-DNDX*DKDY)
      XEST(1,1) = XEST(1,1) - XEST(1,2) 
      YEST(1,1) = YEST(1,1) - YEST(1,2)
      GO TO 30

      END