* Last processed by NICE on 12-Jun-2000 15:54:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
* nic_transform.f
*
      SUBROUTINE TRANSFORM_COORDINATES (U,NU,V,NV,BEAMSEP,ERROR)
*
* fort77 -q +T +FPDVZoui transform.f -L/users/softs/test/lib -lGREG
*
* u       =  input array (r*8) (x,y)
* nu      =  dimension of the input array
* v       =  output array (r*8) (x,y) of dimension
* nv      =  dimension of the output array
* beamsep =  wobbler amplitude (r*8)
*
      INTEGER*4 MAX_SEG, MAX_CROSS, MAX_VTX
      PARAMETER (MAX_VTX = 100, MAX_SEG = 2*MAX_VTX, MAX_CROSS = 50)
      REAL*8 U(MAX_VTX,2), W1(MAX_VTX,2), W2(MAX_VTX,2), V(MAX_SEG,2)
      REAL*8 CRR(MAX_CROSS,2), BOUND(5), GONS(MAX_VTX,4)
      REAL*8 X1, X2, Y1, Y2, X3, Y3, X4, Y4, XC, YC, WS1(MAX_VTX,2)
      REAL*8 W1RAN, W2RAN, W3RAN, W4RAN, WS2(MAX_VTX,2)
      REAL*8 BEAMSEP, W2MIN, XRS, YRS, R, TOT, A, A1, A2, GASDEV
      INTEGER*4 CRP(MAX_CROSS,2), NU, NV, I2MIN, K, L, IS, SA, N
      INTEGER*4 I, J, M, ID, NT, ITER
      LOGICAL IN, GR8_IN, RES, FOUND, ERROR
* Function ATAN2D defined for F2C
C+MAC+LINUX+SOLARIS
      REAL*8  DEG_TO_RAD
      PARAMETER (DEG_TO_RAD=1.74532925199432958D-2)
      REAL*8 XXX,YYY
      ATAN2D(XXX,YYY)=ATAN2(XXX,YYY)/DEG_TO_RAD
C-MAC-LINUX-SOLARIS
      DATA ID /-1/
      SAVE ID
*
      NV = MAX_SEG
      NT = MAX_VTX
      M  = 0
      ITER = 0
      DO WHILE (NV.GT.2*(NU+1)+1.OR.(NV.EQ.NU+1.AND.M.NE.0))
         M = 0
         ITER = ITER+1
         IF (ITER.GT.200) THEN
            ERROR = .TRUE.
            CALL MESSAGE (2,1,'SUPP','Could not transform the support')
            CALL MESSAGE (2,1,'SUPP','You need to redefine a new one')
            RETURN
         ELSE
            ERROR = .FALSE.
         ENDIF
         DO I = 1,NU
            W1RAN = MIN(MAX(GASDEV(ID)*1D-4,-1D-1),1D-1)
            W2RAN = MIN(MAX(GASDEV(ID)*1D-4,-1D-1),1D-1)
            W3RAN = MIN(MAX(GASDEV(ID)*1D-4,-1D-1),1D-1)
            W4RAN = MIN(MAX(GASDEV(ID)*1D-4,-1D-1),1D-1)
            W1(I,1) = U(I,1)+W1RAN+BEAMSEP/2
            W2(I,1) = U(I,1)+W2RAN-BEAMSEP/2
            W1(I,2) = U(I,2)+W3RAN
            W2(I,2) = U(I,2)+W4RAN
         ENDDO
*
* Verify that no point is at the same Y position
*
         FOUND = .TRUE.
         DO WHILE (FOUND)
            FOUND = .FALSE.
            DO I = 1,NU
               IF (W2(I,2).EQ.W1(I,2)) FOUND = .TRUE.
            ENDDO
            IF (FOUND) THEN
               W3RAN = MIN(MAX(GASDEV(ID)*1D-4,-1D-1),1D-1)
               DO I = 1,NU
                  W2(I,2) = U(I,2)+W3RAN
               ENDDO
            ENDIF
         ENDDO
*
* Verify that no point is at the same X position
*
         FOUND = .TRUE.
         DO WHILE (FOUND)
            FOUND = .FALSE.
            DO I = 1,NU-1
               DO J = I,NU
                  IF (W2(I,1).EQ.W1(J,1)) FOUND = .TRUE.
               ENDDO
            ENDDO
            IF (FOUND) THEN
               W2RAN = MIN(MAX(GASDEV(ID)*1D-4,-1D-1),1D-1)
               DO I = 1,NU
                  W2(I,1) = U(I,1)+W2RAN-BEAMSEP/2
               ENDDO
            ENDIF
         ENDDO
*
* Reorganize array ...
*
         I2MIN = 1
         W2MIN = W2(1,1)
         DO I=1,NU
            DO J=1,2
               WS1(I,J) = W1(I,J)
               WS2(I,J) = W2(I,J)
            ENDDO
            IF (W2MIN.NE.W2(I,1)) THEN   ! first point is decisive
               W2MIN = MIN (W2(I,1),W2MIN)
               IF (W2MIN.EQ.W2(I,1)) I2MIN = I
            ENDIF
         ENDDO
         IF (I2MIN.NE.1) THEN
            DO I=1,NU
               K = MOD(I+I2MIN-2,NU)+1
               DO J=1,2
                  W1(I,J) = WS1(K,J)
                  W2(I,J) = WS2(K,J)
               ENDDO
            ENDDO
            I2MIN = 1
            W2MIN = W2(1,1)
         ENDIF
*
* Reset ...
*
         DO J=1,2
            DO I=1,MAX_CROSS
               CRP(I,J) = 0
            ENDDO
         ENDDO
*
* find number of points inside polygons and associated intersections
*
         BOUND(2) =  1D34
         BOUND(3) = -1D34
         BOUND(4) =  1D34
         BOUND(5) = -1D34
         DO I=1,NU
            GONS(I,1) = W2(I,1)
            GONS(I,2) = W2(I,2)
            IF (I.LT.NU) THEN
               GONS(I,3) = W2(I+1,1)-W2(I,1)
               GONS(I,4) = W2(I+1,2)-W2(I,2)
            ELSE
               GONS(I+1,1) = W2(1,1)
               GONS(I+1,2) = W2(1,2)
               GONS(I,3) = W2(1,1)-W2(I,1)
               GONS(I,4) = W2(1,2)-W2(I,2)
            ENDIF
            BOUND(2) = MIN(BOUND(2),W2(I,1)-1D-1)
            BOUND(3) = MAX(BOUND(3),W2(I,1)+1D-1)
            BOUND(4) = MIN(BOUND(4),W2(I,2)-1D-1)
            BOUND(5) = MAX(BOUND(5),W2(I,2)+1D-1)
         ENDDO
         DO I = 1,NU
            IN = GR8_IN (W1(I,1),W1(I,2),NT,NU,GONS,BOUND)
            IF (IN) THEN
               X1 = W1(I,1)
               Y1 = W1(I,2)
               K = MOD(I,NU)+1           ! point after
               X2 = W1(K,1)
               Y2 = W1(K,2)
               DO J = 1,NU
                  L = MOD(J,NU)+1
                  X3 = W2(J,1)
                  Y3 = W2(J,2)
                  X4 = W2(L,1)
                  Y4 = W2(L,2)
                  CALL INTERSEC (X1,Y1,X2,Y2,X3,Y3,X4,Y4,XC,YC,RES)
                  IF (RES) THEN
                     FOUND = .FALSE.
                     DO L = 1,M
                        IF (CRP(L,1).EQ.I+NU.AND.CRP(L,2).EQ.J)
     $                  FOUND = .TRUE.
                     ENDDO
                     IF (.NOT.FOUND) THEN
                        M = M + 1
                        CRR(M,1) = XC
                        CRR(M,2) = YC
                        CRP(M,1) = I+NU
                        CRP(M,2) = J
                     ENDIF
                  ENDIF
               ENDDO
               K = MOD(I+NU-2,NU)+1      ! point before
               X2 = W1(K,1)
               Y2 = W1(K,2)
               DO J = 1,NU
                  L = MOD(J,NU)+1
                  X3 = W2(J,1)
                  Y3 = W2(J,2)
                  X4 = W2(L,1)
                  Y4 = W2(L,2)
                  CALL INTERSEC (X1,Y1,X2,Y2,X3,Y3,X4,Y4,XC,YC,RES)
                  IF (RES) THEN
                     FOUND = .FALSE.
                     DO L = 1,M
                        IF (CRP(L,1).EQ.K+NU.AND.CRP(L,2).EQ.J)
     $                  FOUND = .TRUE.
                     ENDDO
                     IF (.NOT.FOUND) THEN
                        M = M + 1
                        CRR(M,1) = XC
                        CRR(M,2) = YC
                        CRP(M,1) = K+NU
                        CRP(M,2) = J
                     ENDIF
                  ENDIF
               ENDDO
            ENDIF
         ENDDO
*
         BOUND(2) =  1D34
         BOUND(3) = -1D34
         BOUND(4) =  1D34
         BOUND(5) = -1D34
         DO I=1,NU
            GONS(I,1) = W1(I,1)
            GONS(I,2) = W1(I,2)
            IF (I.LT.NU) THEN
               GONS(I,3) = W1(I+1,1)-W1(I,1)
               GONS(I,4) = W1(I+1,2)-W1(I,2)
            ELSE
               GONS(I+1,1) = W2(1,1)
               GONS(I+1,2) = W2(1,2)
               GONS(I,3) = W1(1,1)-W1(I,1)
               GONS(I,4) = W1(1,2)-W1(I,2)
            ENDIF
            BOUND(2) = MIN(BOUND(2),W1(I,1)-1D-1)
            BOUND(3) = MAX(BOUND(3),W1(I,1)+1D-1)
            BOUND(4) = MIN(BOUND(4),W1(I,2)-1D-1)
            BOUND(5) = MAX(BOUND(5),W1(I,2)+1D-1)
         ENDDO
         DO I = 1,NU
            IN = GR8_IN (W2(I,1),W2(I,2),NT,NU,GONS,BOUND)
            IF (IN) THEN
               X1 = W2(I,1)
               Y1 = W2(I,2)
               K = MOD(I,NU)+1           ! point after
               X2 = W2(K,1)
               Y2 = W2(K,2)
               DO J = 1,NU
                  L = MOD(J,NU)+1
                  X3 = W1(J,1)
                  Y3 = W1(J,2)
                  X4 = W1(L,1)
                  Y4 = W1(L,2)
                  CALL INTERSEC (X1,Y1,X2,Y2,X3,Y3,X4,Y4,XC,YC,RES)
                  IF (RES) THEN
                     FOUND = .FALSE.
                     DO L = 1,M
                        IF (CRP(L,1).EQ.J+NU.AND.CRP(L,2).EQ.I)
     $                  FOUND = .TRUE.
                     ENDDO
                     IF (.NOT.FOUND) THEN
                        M = M + 1
                        CRR(M,1) = XC
                        CRR(M,2) = YC
                        CRP(M,1) = J+NU
                        CRP(M,2) = I
                     ENDIF
                  ENDIF
               ENDDO
               K = MOD(I+NU-2,NU)+1      ! point before
               X2 = W2(K,1)
               Y2 = W2(K,2)
               DO J = 1,NU
                  L = MOD(J,NU)+1
                  X3 = W1(J,1)
                  Y3 = W1(J,2)
                  X4 = W1(L,1)
                  Y4 = W1(L,2)
                  CALL INTERSEC (X1,Y1,X2,Y2,X3,Y3,X4,Y4,XC,YC,RES)
                  IF (RES) THEN
                     FOUND = .FALSE.
                     DO L = 1,M
                        IF (CRP(L,1).EQ.J+NU.AND.CRP(L,2).EQ.K)
     $                  FOUND = .TRUE.
                     ENDDO
                     IF (.NOT.FOUND) THEN
                        M = M + 1
                        CRR(M,1) = XC
                        CRR(M,2) = YC
                        CRP(M,1) = J+NU
                        CRP(M,2) = K
                     ENDIF
                  ENDIF
               ENDDO
            ENDIF
         ENDDO
*
*            do i=1,m
*               type *,crp(i,1),crp(i,2),crr(i,1),crr(i,2)
*            enddo
*
* first point, first side
*
*
         V(1,1) = W2MIN
         V(1,2) = W2(I2MIN,2)
         IS = I2MIN                      ! starting side
*
         IF (M.NE.0) THEN
            NV = 1
            DO WHILE ((V(NV,1).NE.V(1,1).OR.V(NV,2).NE.V(1,2)
     $         .OR.IS.EQ.I2MIN).AND.NV.LT.MAX_SEG)
               A1 = 0
               IF (NV.GT.1) THEN
                  XRS = V(NV,1)-V(NV-1,1)
                  YRS = V(NV,2)-V(NV-1,2)
                  IF (YRS.NE.0.OR.XRS.NE.0) A1 = ATAN2D(YRS,XRS)
*                  write (6,'(a,2f10.5,a,f10.5)')
*     _                 'Coord:  ',v(nv,1),v(nv,2),'   Angle: ',a1
               ENDIF
               N = MOD(IS,NU)+1
               TOT = 0
               IF (IS.LE.NU) THEN
                  XRS = W2(N,1)-V(NV,1)
                  YRS = W2(N,2)-V(NV,2)
                  R = XRS**2+YRS**2
                  IF (YRS.NE.0.OR.XRS.NE.0)
     $            A = ABS(ATAN2D(YRS,XRS)-A1)
*                  write (6,'(a,2f10.5,a,f10.5)')
*     _                 'Length: ',xrs,yrs,'   Angle: ',a
                  NV = NV+1
                  SA = -1
                  DO I=1,M
                     IF (CRP(I,2).EQ.IS) THEN
                        XRS = CRR(I,1)-V(NV-1,1)
                        YRS = CRR(I,2)-V(NV-1,2)
                        TOT = XRS**2+YRS**2
                        A2  = A1
                        IF (YRS.NE.0.OR.XRS.NE.0)
     $                  A2  = ABS(ATAN2D(YRS,XRS)-A1)
                        R = MIN(R,TOT)
                        A = MAX(A2,A)
*                        write (6,'(i2,a,i2,a,i2,a,2f10.5,a,f10.5)')
*     _                       i,': (',is,'-',crp(i,1),')',
*     _                       sqrt(r),sqrt(tot),'  Angle: ',a
                        IF (R.EQ.TOT.AND.ABS(A-A2).LT.1E-5) THEN
                           V(NV,1) = CRR(I,1)
                           V(NV,2) = CRR(I,2)
                           SA      = I
                        ENDIF
                     ENDIF
                  ENDDO
                  IF (TOT.EQ.0) THEN
                     V(NV,1) = W2(N,1)
                     V(NV,2) = W2(N,2)
                     IS = IS+1
                  ELSE
                     IF (SA.GE.0) THEN
                        IS = CRP(SA,1)
                        CRP(SA,1) = 0
                        CRP(SA,2) = 0
                     ELSE
                        V(NV,1) = W2(N,1)
                        V(NV,2) = W2(N,2)
                        IS = IS+1
                     ENDIF
                  ENDIF
               ELSE
                  XRS = W1(N,1)-V(NV,1)
                  YRS = W1(N,2)-V(NV,2)
                  R = XRS**2+YRS**2
                  A = A1
                  IF (YRS.NE.0.OR.XRS.NE.0)
     $            A = ABS(ATAN2D(YRS,XRS)-A1)
*                  write (6,'(a,2f10.5,a,f10.5)')
*     _                 'Length: ',xrs,yrs,'   Angle: ',a
                  NV = NV+1
                  SA = -1
                  DO I=1,M
                     IF (CRP(I,1).EQ.IS) THEN
                        XRS = CRR(I,1)-V(NV-1,1)
                        YRS = CRR(I,2)-V(NV-1,2)
                        TOT = XRS**2+YRS**2
                        A2  = A1
                        IF (YRS.NE.0.OR.XRS.NE.0)
     $                  A2 = ABS(ATAN2D(YRS,XRS)-A1)
                        R = MIN(R,TOT)
                        A = MAX(A2,A)
*                        write (6,'(i2,a,i2,a,i2,a,2f10.5,a,f10.5)')
*     _                       i,': (',is,'-',crp(i,2),')',
*     _                       sqrt(r),sqrt(tot),'  Angle: ',a
                        IF (R.EQ.TOT.AND.ABS(A-A2).LT.1E-5) THEN
                           V(NV,1) = CRR(I,1)
                           V(NV,2) = CRR(I,2)
                           SA      = I
                        ENDIF
                     ENDIF
                  ENDDO
                  IF (TOT.EQ.0) THEN
                     V(NV,1) = W1(N,1)
                     V(NV,2) = W1(N,2)
                     IS = IS+1
                  ELSE
                     IF (SA.GE.0) THEN
                        IS = CRP(SA,2)
                        CRP(SA,1) = 0
                        CRP(SA,2) = 0
                     ELSE
                        V(NV,1) = W1(N,1)
                        V(NV,2) = W1(N,2)
                        IS = IS+1
                     ENDIF
                  ENDIF
                  IF (IS.GT.2*NU) IS = IS-NU
               ENDIF
            ENDDO
         ELSE
            DO I=1,NU
               V(I,1) = U(I,1)-BEAMSEP/2
               V(I,2) = U(I,2)
               V(I+NU+1,1) = U(I,1)+BEAMSEP/2
               V(I+NU+1,2) = U(I,2)
            ENDDO
            V(NU+1,1) = V(1,1)
            V(NU+1,2) = V(1,2)
            V(2*NU+2,1) = V(NU+2,1)
            V(2*NU+2,2) = V(NU+2,2)
            NV = 2*(NU+1)
         ENDIF
      ENDDO
      RETURN
      END
*
      SUBROUTINE INTERSEC (X1,Y1,X2,Y2,X3,Y3,X4,Y4,XC,YC,RES)
      REAL*8 WX1, WX2, WY1, WY2, XC, YC, A1, A2, B1, B2
      REAL*8 X1, X2, X3, X4, Y1, Y2, Y3, Y4
      LOGICAL RES
*
      RES = .TRUE.
*
      XC = 0
      YC = 0
*
      WX1 = X1-X2
      WX2 = X3-X4
      WY1 = Y1-Y2
      WY2 = Y3-Y4
      IF (WX1.NE.0.AND.WX2.NE.0) THEN
         A1 = WY1/WX1
         A2 = WY2/WX2
         B1 = Y1-A1*X1
         B2 = Y3-A2*X3
         IF (A1.NE.A2) THEN
            XC = (B2-B1)/(A1-A2)
            YC = A1*XC+B1
            IF (YC.LT.Y1.AND.YC.LT.Y2) RES = .FALSE.
            IF (YC.GT.Y1.AND.YC.GT.Y2) RES = .FALSE.
            IF (YC.LT.Y3.AND.YC.LT.Y4) RES = .FALSE.
            IF (YC.GT.Y3.AND.YC.GT.Y4) RES = .FALSE.
         ELSE
            RES = .FALSE.
         ENDIF
      ELSE
         IF (WX1.EQ.0.AND.WX2.NE.0) THEN
            IF (X1.LT.X3.AND.X1.LT.X4) RES = .FALSE.
            IF (X1.GT.X3.AND.X1.GT.X4) RES = .FALSE.
            IF (RES) THEN
               XC = X1
               A2 = WY2/WX2
               B2 = Y3-A2*X3
               YC = A2*X1+B2
            ENDIF
            IF (YC.LT.Y1.AND.YC.LT.Y2) RES = .FALSE.
            IF (YC.GT.Y1.AND.YC.GT.Y2) RES = .FALSE.
         ENDIF
         IF (WX1.NE.0.AND.WX2.EQ.0) THEN
            IF (X3.LT.X1.AND.X3.LT.X2) RES = .FALSE.
            IF (X3.GT.X1.AND.X3.GT.X2) RES = .FALSE.
            IF (RES) THEN
               XC = X3
               A1 = WY1/WX1
               B1 = Y1-A1*X1
               YC = A1*X3+B1
            ENDIF
            IF (YC.LT.Y3.AND.YC.LT.Y4) RES = .FALSE.
            IF (YC.GT.Y3.AND.YC.GT.Y4) RES = .FALSE.
         ENDIF
         IF (WX1.EQ.0.AND.WX2.EQ.0) RES = .FALSE.
      ENDIF
      RETURN
      END
*
      FUNCTION GASDEV (IDUM)
      INTEGER IDUM, ISET
      REAL*8 FAC, GSET, GASDEV, R, V1, V2, RAN1
      DATA ISET /0/
      SAVE ISET, FAC, GSET, R, V1, V2
      IF (ISET.EQ.0) THEN
         R = 0
         DO WHILE (R.EQ.0.OR.R.GE.1D0)
            V1 = 2*RAN1(IDUM)-1.
            V2 = 2*RAN1(IDUM)-1.
            R = V1**2+V2**2
         ENDDO
         FAC = SQRT(-2*LOG(R)/R)
         GSET = V1*FAC
         GASDEV = V2*FAC
         ISET = 1
      ELSE
         GASDEV = GSET
         ISET = 0
      ENDIF
      RETURN
      END
*
      FUNCTION RAN1 (IDUM)
      INTEGER IX1, IX2, IX3, J
      INTEGER IDUM, M1, M2, M3, IA1, IA2, IA3, IC1, IC2, IC3, IFF
      REAL*8 R(97), RM1, RM2, RAN1
      SAVE IFF, IX1, IX2, IX3, R
      PARAMETER (M1 = 259200, IA1 = 7141, IC1 = 54773, RM1 = 1./M1)
      PARAMETER (M2 = 134456, IA2 = 8121, IC2 = 28411, RM2 = 1./M2)
      PARAMETER (M3 = 243000, IA3 = 4561, IC3 = 51349)
      DATA IFF /0/
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
         IFF = 1
         IX1 = MOD(IC1-IDUM,M1)
         IX1 = MOD(IA1*IX1+IC1,M1)
         IX2 = MOD(IX1,M2)
         IX1 = MOD(IA1*IX1+IC1,M1)
         IX3 = MOD(IX1,M3)
         DO J=1,97
            IX1 = MOD(IA1*IX1+IC1,M1)
            IX2 = MOD(IA2*IX2+IC2,M2)
            R(J) = (FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
         ENDDO
         IDUM = 1
      ENDIF
      IX1 = MOD(IA1*IX1+IC1,M1)
      IX2 = MOD(IA2*IX2+IC2,M2)
      IX3 = MOD(IA3*IX3+IC3,M3)
      J = 1+(97*IX3)/M3
      IF (J.GT.97.OR.J.LT.1) J = 97
      RAN1 = R(J)
      R(J) = (FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
      END