C     Last change:  PG   20 Feb 2002
      SUBROUTINE CalcShad(task,
     &     sol_dir_deg,sol_elev_deg,
     &     Coord_Surface,
     &     N_Shades, N_Points, Coord_Shades, Transp_Shades,
     &     Transp_Result,Three,Four)

* Updates
* 020220 PG Revision of area & transparency caclulation
*           Simplified calc. projection matrix
* 980315 AB Sun's position [deg & rad], in separate variables
* 980306 ML Major revision

C
      DOUBLE PRECISION pi
      PARAMETER(pi = 3.141592653589793)
C
cabg  DOUBLE PRECISION sol_dir, sol_elev
      DOUBLE PRECISION sol_dir_deg, sol_elev_deg
      DOUBLE PRECISION Coord_Surface(3,4)
      INTEGER          task, N_Shades, N_points, Three, Four
      DOUBLE PRECISION Coord_Shades(3,N_Shades*4)
      DOUBLE PRECISION Transp_Shades(N_Shades)
      DOUBLE PRECISION Transp_Result
C
C * Inputs
C     task           - 1 for direct,
C                      2 for diffuse, which allows sun below horizon
C     sol_dir_deg    - Direction to the sun [deg]
C     sol_elev_deg   - Elevation of the sun [deg]
C     Coord_Surface  - Coordinates of the surface to project upon
C     N_Shades       - Number of shades (rectangles)
C     Coord_Shades   - Coordinates for the shades
C     Transp_Shades  - Transparency of each shade [0.0, 1.0]
C                        1.0 = transparent
C * Outputs:
C     Transp_Result  - Result as equivalent transparency [0.0,1.0]
C                        1.0 = transparent
C
      INTEGER I, Ii, Iii, nn_shades
      INTEGER, DIMENSION(N_Shades) :: shade_index
      DOUBLE PRECISION sol_dir, sol_elev
      DOUBLE PRECISION, DIMENSION(4,N_Shades*4) ::  Coord_Proj
      DOUBLE PRECISION, DIMENSION(4,4) ::  Coord_Surf_Proj
      DOUBLE PRECISION Trans(4,4)
      DOUBLE PRECISION Translat(4,4)
C-    DOUBLE PRECISION X_Rot, Y_rot, Z_rot
C ??  DOUBLE PRECISION Sol_Vec(4), R, Normal_to_surf(4)
      DOUBLE PRECISION Sol_Vec(4), R, Rx, Ry, Normal_to_surf(4)
C=    DOUBLE PRECISION Surf_X(4), Surf_Y(4), Surf_tmp(4)
      DOUBLE PRECISION Surf_X(4), Surf_Y(4)
      DOUBLE PRECISION Scalar_Product
      DOUBLE PRECISION C1(3),C2(3),C3(3),C4(3)
      LOGICAL Debug,debug1

      DOUBLE PRECISION P0(4),Px(4),Py(4), Rotat(4,4)

C
C * Initializations
C
      Debug = .FALSE.
      Debug1 = .FALSE.
      IF(DEBUG) THEN
        WRITE(*,*) '==============='
        Write(*,*) 'S ',sol_dir_deg, sol_elev_deg
      END IF

      IF(DEBUG1)
     & Write(127,'(a,2f10.2,a)') '''(',sol_dir_deg, sol_elev_deg,')'
C
C * Transform Sun position to a vector
C
      Sol_dir = sol_dir_deg * PI / 180.0
      sol_elev = sol_elev_deg * PI / 180.0
      Sol_Vec =
     &  (/ Sin(Sol_Dir) * Cos(Sol_Elev),
CCCC switch!
     &     Cos(Sol_Dir) * Cos(Sol_Elev),
     &     Sin(Sol_Elev),
     &     1.0D0 /)
c>      Debug = .TRUE.
      Transp_Result = 0.0D0

C * If the sun is under horizon - return unless calc diffuse shading
cle>  IF(sol_elev_deg .LE. 0.0D0 .AND. task .EQ. 0) RETURN
      IF(sol_elev_deg .LE. 0.0D0 .AND. task .EQ. 1) RETURN
C
C 1 3, 4 3
C It is only interesting to align window with XY, hence the Abs
C ?? Very strange order of points in Coord_Surf:
C ??  4 2
C ??  3 1
C ??  Find an ortogonal coordinate system where
C ??  3 --> (0,0,0)
C ??  1 --> (256, 0, 0)
C ??  4 --> (0, 256, 0)
C ??  and the transformation matrix to such coordinates
C 
C
      Surf_X(1:3) = Abs(Coord_Surface(:,1) - Coord_Surface(:,3))
C ??  Surf_X(1:3) = Coord_Surface(:,1) - Coord_Surface(:,3)
      Surf_X(4) = 1.0
      Surf_Y(1:3) = Abs(Coord_Surface(:,4) - Coord_Surface(:,3))
C ??  Surf_Y(1:3) = Coord_Surface(:,4) - Coord_Surface(:,3)
      Surf_Y(4) = 1.0
C
      CALL Cross_product(Normal_To_Surf, Surf_Y, Surf_X)
C ??  CALL Cross_product(Normal_To_Surf, Surf_X, Surf_Y)
C fix
C ?? this is wrong if the vector has coordinates with different signs ??
      Normal_to_surf = ABS(Normal_to_surf)
      IF(DEBUG) THEN
        Write(*,'(1x,A,1x,3F6.2)') 'Sol_Vec',Sol_Vec(1:3)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,1)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,2)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,3)
        Write(*,'(1x,A,1x,3F6.2)') 'Coord_surface',Coord_surface(1:3,4)
        Write(*,'(1x,A,1x,3F6.2)') 'Normal_To_Surf',Normal_To_Surf(1:3)
        Write(*,'(1x,A,1x,3F6.2)') 'Scalar_product',
     & Scalar_product(Sol_Vec, Normal_To_Surf)
      ENDIF
      IF (Scalar_product(Sol_Vec, Normal_To_Surf) .LE. 0.0D0) RETURN



C
C * Make transformation matrix for the surface on the XY-plane
C
C     * Normalize the surface x and y
C
      Rx = Sqrt(
     &     Surf_X(1)*Surf_X(1) +
     &     Surf_X(2)*Surf_X(2) +
     &     Surf_X(3)*Surf_X(3))
      Surf_X(1:3) = Surf_X(1:3)/Rx
      Ry = Sqrt(
     &     Surf_Y(1)*Surf_Y(1) +
     &     Surf_Y(2)*Surf_Y(2) +
     &     Surf_Y(3)*Surf_Y(3))
      Surf_Y(1:3) = Surf_Y(1:3)/Ry
C
C     * Transformation matrix
C
C     * Tranclation
C
      Translat = 0.0
      DO i = 1,4
         Translat(i,i) = 1.0
      END DO

      Translat(1:3,4) = - Coord_Surface(:,3)

C
C     * Rotation (Surf_X and Surf_Y should be ortogonal!)
C
      Rotat = 0
      Rotat(1,1:3) = Surf_X(1:3)
      Rotat(2,1:3) = Surf_Y(1:3)
      Rotat(3,1:3) = Normal_To_Surf(1:3)
C      Rotat(4,1:3) = 0

      Trans = MATMUL(Rotat,Translat)


      P0(1:3) = Coord_Surface(1:3,3)
      Px(1:3) = Coord_Surface(1:3,1)
      Py(1:3) = Coord_Surface(1:3,4)
      P0(4) = 1
      Px(4) = 1
      Py(4) = 1

      IF(Debug)THEN
         Write(*,'(A,/4(4f6.2/)/)') ' Translat ',(Translat(:,i),i=1,4)
         Write(*,'(A,/4(4f6.2/)/)') ' Rotat   ',(Rotat(:,i),i=1,4)
         Write(*,'(A,/4(4f6.2/)/)') ' Trans    ',(Trans(:,i),i=1,4)
         Write(*,'(A,4f6.2)') ' trans on 0 ',MATMUL(Trans, P0)
         Write(*,'(A,4f6.2)') ' trans on X ',MATMUL(Trans, Px)
         Write(*,'(A,4f6.2)') ' trans on Y ',MATMUL(Trans, Py)
      END IF
C      CALL NMF_ERROR("Exit")
C
C     * Transform the surface
C
      Coord_Surf_Proj(1:3,:) = Coord_Surface
      Coord_Surf_Proj(4,:) = 1.0
      Coord_Surf_Proj = MATMUL(Trans, Coord_Surf_Proj)
      IF(Debug)THEN
         Write(*,*) 'Projected surface'
         Write(*,*) Nint(Coord_Surf_Proj(1:2,1)*100)/100.0
         Write(*,*) Nint(Coord_Surf_Proj(1:2,2)*100)/100.0
         Write(*,*) Nint(Coord_Surf_Proj(1:2,3)*100)/100.0
         Write(*,*) Nint(Coord_Surf_Proj(1:2,4)*100)/100.0
      END IF
C
C * Transform every Shade by transforming the coordinate vector
C
      Coord_Proj(1,:) = Coord_Shades(1,:)
      Coord_Proj(2,:) = Coord_Shades(2,:)
      Coord_Proj(3,:) = Coord_Shades(3,:)
      Coord_Proj(4,:) = 1.0
      IF(Debug)THEN
        write(*,*) 'Coord_proj before trans'
        do i=1,n_shades*4
          Write(*,*) Coord_proj(1:3,I)
        end do
      END IF
C
      Coord_Proj = MATMUL(Trans, Coord_Proj)
C
      IF(Debug)THEN
        write(*,*) 'Coord_proj after trans'
        do i=1,n_shades*4
          Write(*,*) Coord_proj(1:3,I)
        end do
      END IF
C
C * Cut with XY-plane
C
C ?? The intersection may be pentagonal ??

      nn_shades = 0
      DO I=1,n_shades
        Ii = (I-1)*4
        C1=Coord_proj(1:3,Ii+1)
        C2=Coord_proj(1:3,Ii+2)
        C3=Coord_proj(1:3,Ii+3)
        C4=Coord_proj(1:3,Ii+4)
Cabg?? Ok to update on the fly?
        call CutZZero(C1,C2)
        call CutZZero(C2,C3)
        call CutZZero(C3,C4)
        call CutZZero(C4,C1)
        Coord_proj(1:3,Ii+1)=C1
        Coord_proj(1:3,Ii+2)=C2
        Coord_proj(1:3,Ii+3)=C3
        Coord_proj(1:3,Ii+4)=C4
        IF(C1(3) .LE. 0 .AND. C2(3) .LE. 0 .AND.
     &     C3(3) .LE. 0 .AND. C4(3) .LE. 0) THEN
C * Discard
        ELSE
C * Cut
           nn_shades=nn_shades+1 
           shade_index(nn_shades)=I
        END IF
      end do

      do i=1,nn_shades
        Ii = (i-1)*4
        Iii = (shade_index(i)-1)*4
        Coord_proj(1:3,Ii+1)=Coord_proj(1:3,Iii+1)
        Coord_proj(1:3,Ii+2)=Coord_proj(1:3,Iii+2)
        Coord_proj(1:3,Ii+3)=Coord_proj(1:3,Iii+3)
        Coord_proj(1:3,Ii+4)=Coord_proj(1:3,Iii+4)
      end do
      n_shades = nn_shades
C
      IF(Debug)THEN
        write(*,*) 'Coord_proj efter CutZZero'
        do i=1,n_shades*4
          Write(*,*) Coord_proj(1:3,I)
        end do
      END IF
C
C * Only rotation of the solar direction
C
      Sol_vec(1:3) = MATMUL(Trans(1:3,1:3),Sol_vec(1:3))
      IF(Debug)THEN
C-       Write(*,*) 'X_Rot,Y_rot, Sol_Dir, Sol_Elev'
C-       Write(*,*) X_Rot,Y_rot, Sol_Dir, Sol_Elev
         Write(*,*) 'Sol_Vec'
         Write(*,'(4f6.2)') Sol_Vec
      END IF
C
C
C * Project every point in the shades to the transformed surface (Z=0)
C
      IF(Debug)THEN
         Write(*,*) 'Coord_Proj'
         IF(debug1) write(127,*) '(mapcar #''draw-1 ''('
      END IF
      DO I = 1, 4*N_Shades
         R=Coord_Proj(3,I)/Sol_Vec(3)
         Coord_Proj(1:3,I) = Coord_Proj(1:3,I) - R * Sol_Vec(1:3)
         Coord_Proj(1:3,I) = Nint(Coord_Proj(1:3,I)*100) / 100.0
         IF(Debug)THEN
            Write(*,'(1x,2F7.2)') Coord_Proj(1:2,I)
         END IF
         IF(Debug1)THEN
            if(mod(I,4).eq.1) write(127,*) ' ('
            write(127,'(3x,a1,3F6.2,a1)') '(',Coord_Proj(1:3,I),')'
            if(mod(I,4).eq.0) write(127,*) ' )'
         END IF
      END DO
      IF(Debug1)THEN
         write(127,*) '))'
      END IF
C
C * Scanning
C
      CALL Scan_and_cut(
     &     Coord_surf_proj, N_Shades, Coord_proj,
     &     Transp_Shades, Transp_Result)
C>      Debug = .TRUE.
      IF(DEBUG) WRITE(*,*) Transp_Result
      RETURN
      END



***************************************************************
*
*	  Calculate the transparenty of a system of 2d-shades
*
C
C * Inputs
C     Coord_Surface    - Coordinates of the surface to project upon
C     N_Shades         - Number of shades (rectangles)
C     Coord_Shades     - Coordinates for the shades
C     Transp_Shades    - Transparency of each shade [0.0, 1.0]
C                        1.0 = transparent
C * Outputs:
C     Transp_Result    - Result as equivalent transparency [0.0,1.0]
C


      SUBROUTINE Scan_and_cut(
     &     Coord_surface, N_Shades, Coord_shades,
     &     Transp_Shades, Transp_Result)
C
      DOUBLE PRECISION Coord_Surface(4,4)
      INTEGER N_Shades
      DOUBLE PRECISION Coord_Shades(4,N_Shades*4)
      DOUBLE PRECISION Transp_Shades(N_Shades)
      DOUBLE PRECISION Transp_Result
      LOGICAL Debug
      INTEGER Ip1, Ip2, Ip3, Ip4
      DOUBLE PRECISION Xmin, Xmax, Ymin, Ymax, Total_area
      DOUBLE PRECISION Xp1, Xp2, Xp3, Xp4, Yp1, Yp2, Yp3, Yp4
      DOUBLE PRECISION Eps
      PARAMETER (Eps=1.0D-8)
C
C * Pointers etc.
C
      INTEGER Nullp, SUCC, PRE, DAT, TRANS, AREA
      PARAMETER (Nullp=-1,
     &           SUCC=1, PRE=2, DAT=3, TRANS=1, AREA=2)
      INTEGER MAX_Surf, MAX_Pnt, Max_Coord
      PARAMETER (MAX_Surf=1000,
     &           MAX_Pnt=4000,
     &           Max_Coord=4000)
      INTEGER M_Surf, M_Pnt, M_Coord
      INTEGER Slistp(Max_Surf), Slist(DAT,Max_Surf)
      INTEGER Clistp(Max_Coord)
      DOUBLE PRECISION Clist(2,Max_Coord)
      DOUBLE PRECISION RSlist(AREA,Max_Surf)
      INTEGER Plistp(Max_Pnt), Plist(3,Max_Pnt)
C Functions
      INTEGER New_Coord, New_Surf
C-    DOUBLE PRECISION Trans_surfaces
C
C=    INTEGER Root_out, Rootp, Root_scan
      INTEGER Rootp
C=    INTEGER I_SurfP, I_Coord, I_Surfout, I_Scan
      INTEGER I_SurfP, I_Coord
      INTEGER I_P1, I_P2, I_PS, I, I_PX
      DOUBLE PRECISION Coord(2), Coord1(2), Coord2(2), Coordx(2), R
      DOUBLE PRECISION CoordQ1(2), CoordQ2(2)
C=    INTEGER I_SurfQ, I_Q1, I_Q2, I_PScan
      INTEGER I_SurfQ, I_Q1, I_Q2
      LOGICAL intersect, Lmin
C-    DOUBLE PRECISION Trans_w_size

C+   [PG 20/02/2002]
      INTEGER N_X, N_Shad, N_Edges
	  INTEGER relpos, j, k, jj, js
	  DOUBLE PRECISION xLeft, xRight, X, X1, X2, Y, Y0, TRAN, YTRAN
	  DOUBLE PRECISION XARR(Max_Coord), YARR(2*MAX_Surf)
	  DOUBLE PRECISION TRARR(MAX_Surf)
	  INTEGER XORD(Max_Coord), YORD(2*MAX_Surf), SHARR(2*MAX_Surf)
	  LOGICAL ONARR(MAX_Surf)
C.   [PG 20/02/2002]

C
C * SList (1,:) = SUCC,  - list of shades
C *       (2,:) = PRE,
C *       (3,:) = DAT -> PList,  NullP = -1
C
C * PList (1,:) = SUCC,	 - list of points (in shade)
C *       (2,:) = PRE,
C *       (3,:) = DAT -> CList,  NullP = -1
C
C * CList (1,:) = X
C * CList (2,:) = Y
C
C *RSList (1,:) = TRANSP,
C *       (2,:) = AREA,
C
C * Initializations
C
      Debug = .FALSE.
      IF (DEBUG) THEN
          WRITE(*,*) 'Scan_and_cut'
      ENDIF
C
      Rootp = Nullp
C-    Root_out = Nullp
C-    Root_scan = Nullp
C      
      DO I = 1, Max_Surf
          Slistp(I)=I
      END DO
      M_Surf=0

      DO I = 1, Max_Coord
          Clistp(I)=I
      END DO
      M_Coord=0

      DO I = 1, Max_Pnt
          Plistp(I)=I
      END DO
      M_Pnt=0
C
C-    Trans_w_size = 0.0D0
C ======================================================================
C
C * Move every shade that is potentially inside the surface to
C * a list of shadeareas SList
C
      Xmin = Coord_Surface(1,1)
      Xmax = Coord_Surface(1,3)
      Ymin = Coord_Surface(2,1)
      Ymax = Coord_Surface(2,3)
      DO I=1,4
         Xmin = Min(Xmin,Coord_Surface(1,I))
         Xmax = Max(Xmax,Coord_Surface(1,I))
         Ymin = Min(Ymin,Coord_Surface(2,I))
         Ymax = Max(Ymax,Coord_Surface(2,I))
      END DO
      Total_area = (Xmax - Xmin) * (Ymax - Ymin)
      DO I=1,N_Shades
        Ip1 = (I-1)*4+1
        Ip2 = (I-1)*4+2
        Ip3 = (I-1)*4+3
        Ip4 = (I-1)*4+4
        Xp1 = Coord_Shades(1,Ip1)
        Yp1 = Coord_Shades(2,Ip1)
        Xp2 = Coord_Shades(1,Ip2)
        Yp2 = Coord_Shades(2,Ip2)
        Xp3 = Coord_Shades(1,Ip3)
        Yp3 = Coord_Shades(2,Ip3)
        Xp4 = Coord_Shades(1,Ip4)
        Yp4 = Coord_Shades(2,Ip4)
        IF (
     &     (Xmin .GT. XP1 .AND. Xmin .GT. XP2 .AND.
     &      Xmin .GT. XP3 .AND. Xmin .GT. XP4)) THEN
C Discard left
CC           Write(*,*) 'l',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE IF (
     &     (Xmax .LT. XP1 .AND. Xmax .LT. XP2 .AND.
     &      Xmax .LT. XP3 .AND. Xmax .LT. XP4)) THEN
C Discard right
CC           Write(*,*) 'r',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE IF (
     &     (Ymin .GT. YP1 .AND. Ymin .GT. YP2 .AND.
     &      Ymin .GT. YP3 .AND. Ymin .GT. YP4)) THEN
C Discard below
CC           Write(*,*) 'b',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE IF (
     &     (Ymax .LT. YP1 .AND. Ymax .LT. YP2 .AND.
     &      Ymax .LT. YP3 .AND. Ymax .LT. YP4)) THEN
C Discard above
CC           Write(*,*) 'a',Xp1,yp1
CC           Write(*,*) Xp2,yp2
CC           Write(*,*) Xp3,yp3
CC           Write(*,*) Xp4,yp4
        ELSE
C Do not discard!
          I_SurfP = New_surf(M_Surf,Rootp,Slistp,Slist,RSlist,
     &          Transp_Shades(I))
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip1))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip2))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip3))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          I_Coord = New_coord(M_Coord,Clistp,Clist,Coord_Shades(1,Ip4))
          Call Add_coord(I_SurfP,I_Coord,Slist,M_Pnt,Plistp,Plist)
          IF (debug) THEN
            WRITE(*,*) 'Added shade'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF
      END DO
C
C * Now we have all surfaces stored in Slist beginning from Rootp
C
C ======================================================================
C
C * Start by cutting against the window surface
C * 1st left & right against Xmin & Xmax
C
      I_SurfP = Rootp
 2000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 2050
      I_P1=Slist(DAT,I_SurfP)
C>ML:980904 moved down>      Coord1(:) = Clist(:,Plist(DAT,I_P1))
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 2000
      END IF
C>>
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 2100 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
C cut Left and Right
      IF( Abs(Coord1(1)-Coord2(1)) .GT. Eps ) THEN
        Lmin = Coord1(1) .LE. Coord2(1)
        I_PX = I_P1
CC      write(*,*) 'xm',Xmin,Xmax
CC      write(*,*) 'coord1',coord1
CC      write(*,*) 'coord2',coord2
        R=(Xmin-Coord1(1))/(coord2(1)-Coord1(1))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(1)=Xmin+eps
          Coordx(2)=R*(Coord2(2)-Coord1(2))+Coord1(2)
          IF(Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)-(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_PX=Plist(SUCC,Plist(SUCC,I_P1))
          ELSE
             Coord(:)=Coordx(:)-(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '2100 cut left'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
        R=(Xmax-Coord1(1))/(coord2(1)-Coord1(1))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(1)=Xmax-eps
          Coordx(2)=R*(Coord2(2)-Coord1(2))+Coord1(2)
          IF(.NOT. Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)+(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          ELSE
             Coord(:)=Coordx(:)+(/2*eps,0.0D0/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF       
          IF (debug) THEN
            WRITE(*,*) '2100 cut right'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF
CC        write(*,*) 'coords lr',coordx
      END IF
      Coord1(:)=Coord2(:)
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 2100
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 2000
 2050 CONTINUE
C>
CC      CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C ----------------------------------------------------------
C Remove right and left points
      I_SurfP = Rootp
 3000 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 3050
      I_P1=Slist(DAT,I_SurfP)
 3100 CONTINUE
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 3000
      END IF
      I_Coord = Plist(DAT,I_P1)
      Coord1(:) = Clist(:,I_Coord)
C cut Left and Right
      IF( Coord1(1) .LT. Xmin ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3100
        END IF
      ELSE IF( Coord1(1) .GT. Xmax ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3100
        END IF
      ELSE
      END IF
      I_P1 = Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) GOTO 3100
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 3000
 3050 CONTINUE
      IF (debug) THEN
        WRITE(*,*) '3050 removed left & right'
        CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
      END IF
C ======================================================================
C
C * Now bottom & top against Ymin & Ymax
C
      I_SurfP = Rootp
 2200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 2250
      I_P1=Slist(DAT,I_SurfP)
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 2200
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 2300 CONTINUE
      I_P2=Plist(SUCC,I_P1)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
C cut Up and down
      IF( Abs(Coord1(2)-Coord2(2)) .GT. Eps ) THEN
        Lmin = Coord1(2) .LE. Coord2(2)
        I_PX = I_P1
C      write(*,*) 'ym',Ymin,Ymax
C      write(*,*) 'coord1',coord1
C      write(*,*) 'coord2',coord2
        R=(Ymin-Coord1(2))/(coord2(2)-Coord1(2))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(2)=Ymin+eps
          Coordx(1)=R*(Coord2(1)-Coord1(1))+Coord1(1)
          IF(Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)-(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_PX=Plist(SUCC,Plist(SUCC,I_P1))
          ELSE
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)-(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '2300 cut bottom'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
        R=(Ymax-Coord1(2))/(coord2(2)-Coord1(2))
        IF (R .GE. 0.0 .AND. R+eps .LT. 1.0) THEN
          Coordx(2)=Ymax-eps
          Coordx(1)=R*(Coord2(1)-Coord1(1))+Coord1(1)
          IF(.NOT. Lmin)THEN
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             Coord(:)=Coordx(:)+(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          ELSE
             Coord(:)=Coordx(:)+(/0.0D0,2*eps/)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
             I_Coord = New_coord(M_Coord,Clistp,Clist,Coordx)
             Call Add_coord_after(I_PX,I_Coord,Slist,M_Pnt,Plistp,Plist)
          END IF
          IF (debug) THEN
            WRITE(*,*) '2300 cut top'
            CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
          END IF
        END IF       
CC        write(*,*) 'coords bt',coordx
      END IF
C      
      Coord1(:)=Coord2(:)
      I_P1 = I_P2
      IF (I_P1 .NE. Nullp) GOTO 2300
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 2200
 2250 CONTINUE
C>
C      CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C ----------------------------------------------------------
C
C Remove top and bottom points
C
      I_SurfP = Rootp
 3200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 3250
      I_P1=Slist(DAT,I_SurfP)
 3300 CONTINUE
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 3200
      END IF
      I_Coord = Plist(DAT,I_P1)
      Coord1(:) = Clist(:,I_Coord)
C cut Top and Down
      IF( Coord1(2) .LT. Ymin ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3300
        END IF
      ELSE IF( Coord1(2) .GT. Ymax ) THEN
C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
        I_PS = Plist(PRE,I_P1)
        CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
        IF( I_PS .NE. Nullp)THEN
          I_P1 = I_PS
        ELSE
          I_P1=Slist(DAT,I_SurfP)
          GOTO 3300
        END IF
      ELSE
      END IF
      I_P1 = Plist(SUCC,I_P1)
      IF (I_P1 .NE. Nullp) GOTO 3300
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 3200
 3250 CONTINUE
C
C>
      IF(Debug)THEN
        write(*,*) 'Rootp scan'
        CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C-        CALL Remove_duplicates (Rootp,
C-   &  Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
      END IF
C
C ======================================================================
C
C * Remove void surfaces
C * Add intersects between surfaces
C   [PG 20/02/2002] Make an array of x-coordinates
C                   The intersections are added to this array
C

C+ [PG 20/02/2002]
      N_X = 0
C. [PG 20/02/2002]

      I_SurfP = Rootp
 7400 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 7450
      I_P1=Slist(DAT,I_SurfP)
      IF(I_P1 .EQ. Nullp) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 7400
      END IF
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
 7500 CONTINUE
C+ [PG 20/02/2002]
	  N_X = N_X+1
	  XARR(N_X) = Coord1(1)
C. [PG 20/02/2002]
      I_P2=Plist(SUCC,I_P1)
cle>  Coord2(:) = Clist(:,Plist(DAT,I_P2))
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF

      I_SurfQ = Slist(SUCC,I_SurfP)
 7600 CONTINUE
      IF(I_SurfQ .EQ. Nullp) GOTO 7650
      I_Q1=Slist(DAT,I_SurfQ)
      IF (I_Q1 .EQ. Nullp) THEN
        CALL Remove_surf(I_SurfQ,Rootp,M_Surf,Slistp,Slist)
        I_SurfQ = Slist(SUCC,I_SurfQ)
        GOTO 7600
      END IF
      CoordQ1(:) = Clist(:,Plist(DAT,I_Q1))
 7700 CONTINUE
      I_Q2=Plist(SUCC,I_Q1)
      IF (I_Q2 .EQ. Nullp) THEN
        CoordQ2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfQ)))
      ELSE
        CoordQ2(:) = Clist(:,Plist(DAT,I_Q2))
      ENDIF
      IF (intersect(Coord1,Coord2,CoordQ1,CoordQ2,Coord)) THEN
C- [PG 20/02/2002]
C-         I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
C-         Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-         I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
C-         Call Add_coord_after(I_Q1,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-         IF (debug) THEN
C-           WRITE(*,*) '7700 intersect surfaces'
C-           CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C-         END IF
C+ [PG 20/02/2002]
	  N_X = N_X+1
	  XARR(N_X) = Coord(1)
C. [PG 20/02/2002]
      END IF
      IF (I_Q2 .NE. Nullp) THEN
        I_Q1=I_Q2
        CoordQ1 = CoordQ2
        GOTO 7700
      END IF
      I_SurfQ = Slist(SUCC,I_SurfQ)
      GOTO 7600
 7650 CONTINUE
C
      IF (I_P2 .NE. Nullp) THEN
        I_P1=I_P2
        Coord1 = Coord2
        GOTO 7500
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7400
 7450 CONTINUE
C
C ----------------------------------------------------------
C [PG 20/02/20002] New scan algorithm
C 

C Initialize
      Transp_result = 0.0

C * Sort the array of x-coordinates

      N_X = N_X+1
	  XARR(N_X) = Xmax
      CALL ORDER(XARR,XORD,N_X)

C * Scan stripes

      xRight = Xmin
      DO 2290 i=1,N_X
      xLeft = xRight
	  xRight = XARR(XORD(I))
	  IF (xRight-xLeft<Eps) GOTO 2290
	  X = 0.5*(xLeft+xRight)

C * Scan all shades at level x

      N_Shad = 0
	  N_Edges = 0
      I_SurfP = Rootp
 2220 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 2225
      I_P1=Slist(DAT,I_SurfP)
 2222 CONTINUE
      IF(I_P1 .EQ. Nullp .OR. RSList(1,I_SurfP) .LT. 0) THEN
         I_SurfP = Slist(SUCC,I_SurfP)
         GOTO 2220
      END IF

C Start scanning shade
      Coord1(:) = Clist(:,Plist(DAT,I_P1))
	  X1 = Coord1(1)
      I_P2 = I_P1

C mark position: -1 to left of x; 0 - intersects; 1 - to right
      IF (X1 .GT. X) THEN
         relpos = 1
	  ELSE
	     relpos = -1
	  END IF

 2221 CONTINUE
      I_P2 = Plist(SUCC,I_P2)
      IF (I_P2 .EQ. Nullp) THEN
         Coord2(:) = Clist(:,Plist(DAT,I_P1))
      ELSE
         Coord2(:) = Clist(:,Plist(DAT,I_P2))
      ENDIF
	  X2 = Coord2(1)
C  Skip vertical lines
	  IF(ABS(X1-X2) .LE. Eps ) GOTO 2224
C  Find sides that intersects X	 
      R = (X-X1)/(X2-X1)
      IF (R .LT. 0.0 .OR. R .GT. 1.0) GOTO 2224
C  For first intersection: store dhade data
	  IF (relpos .NE. 0) THEN
	    relpos = 0
		N_Shad = N_Shad + 1
		TRARR(N_SHAD) = RSList(1,I_SurfP)
  	    ONARR(N_SHAD) = .FALSE.
	  END IF
C  Store intersection point
      N_Edges = N_Edges+1
	  YARR(N_Edges) = R*(Coord2(2)-Coord1(2))+Coord1(2)
	  SHARR(N_Edges) = N_Shad
C  Next edge or shade
 2224 CONTINUE
	  IF (I_P2 .EQ. Nullp) THEN
C    Shades to left of x: do not scan again
	    IF (relpos .LT. 0) THEN
		   RSList(1,I_SurfP) = -1
		END IF
        I_SurfP = Slist(SUCC,I_SurfP)
        GOTO 2220
	  END IF
	  Coord1 = Coord2
	  X1 = X2
	  GOTO 2221

 2225 CONTINUE
C  Calculate equivalent transparency in this stripe
      Y0 = Ymin
	  TRAN = 1.0
	  YTRAN = 0.0
      IF (N_Shad .EQ. 0) GOTO 2235
C  Scan y's
	  CALL ORDER(YARR, YORD, N_Edges)
	  DO 2230 j=1,N_Edges
	  jj =  YORD(j)
	  Y = YARR(jj)
	  js = SHARR(jj)
	  YTRAN = YTRAN + (Y-Y0)*TRAN
	  Y0 = Y
	  IF (ONARR(js)) THEN
C   end of shade: switch off and recalculate the transparency
        ONARR(js) = .FALSE.
		TRAN = 1
		DO 2231 k=1,N_Shad
		  IF (ONARR(k)) TRAN = TRAN*TRARR(k)
 2231   CONTINUE
      ELSE
C   start of shade: switch on and multiply the transparency
        ONARR(js) = .TRUE.
		TRAN = TRAN*TRARR(js)
	  END IF
 2230 CONTINUE

 2235 CONTINUE
      YTRAN = YTRAN + (YMax-Y0)

C  Add to overall transparency

      Transp_result = Transp_result + (xRight-xLeft)*YTRAN

C      WRITE	(*,*) 'Strip: ', xLeft, ' - ', xRight
C	  WRITE (*,*) 'y: ', N_Edges, (YARR(k), k=1,N_Edges)
C	  WRITE (*,*) 'Shades: ', N_Shad, (SHARR(k), k=1,N_Shad)
C	  WRITE (*,*) 'Transp: ', YTRAN, Transp_result

C  Next stripe
 2290 CONTINUE

C  Result

      Transp_result = Transp_result / Total_area
	  RETURN
	  END
      

C- 
C- C
C- C ----------------------------------------------------------
C- C
C- C * Combine all surfaces to one scan list
C- C
C-       I_SurfP = Rootp
C-       I_Scan = New_surf(M_Surf,Root_scan,Slistp,Slist,RSlist,0.0D0)
C-  7000 CONTINUE
C-       IF(I_SurfP .EQ. Nullp) GOTO 7050
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF(I_P1 .EQ. Nullp) THEN
C-          I_SurfP = Slist(SUCC,I_SurfP)
C-          IF(I_SurfP .EQ. Nullp) GOTO 7050
C- C> ????
C- C>         GOTO 7400
C-          GOTO 7000
C-       END IF
C-       Coord1(:) = Clist(:,Plist(DAT,I_P1))
C-       I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
C-       Call Add_coord(I_Scan,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-       IF (debug) THEN
C-         WRITE(*,*) '7000 combine'
C-         CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
C-       END IF
C-  7100 CONTINUE
C-       I_P1=Plist(SUCC,I_P1)
C-       IF (I_P1 .NE. Nullp) THEN
C-         Coord1(:) = Clist(:,Plist(DAT,I_P1))
C-         I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
C-         Call Add_coord(I_Scan,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-         IF (debug) THEN
C-           WRITE(*,*) '7100 combine'
C-           CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
C-         END IF
C-         GOTO 7100
C-       END IF
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 7000
C-  7050 CONTINUE
C- 
C-       IF (debug) THEN
C-         WRITE(*,*) 'before sort'
C-         CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
C-       END IF
C- C
C- C ======================================================================
C- C
C- C Sort points and remove duplicate x-values
C- C
C-  7900 CONTINUE
C-       I_SurfP = Root_scan
C-  8000 CONTINUE
C-       IF(I_SurfP .EQ. Nullp) GOTO 8050
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF(I_P1 .EQ. Nullp) THEN
C-          I_SurfP = Slist(SUCC,I_SurfP)
C-          GOTO 8000
C-       END IF
C-       Coord1(:) = Clist(:,Plist(DAT,I_P1))
C-  8100 CONTINUE
C-       I_P2=Plist(SUCC,I_P1)
C-       IF (I_P2 .NE. Nullp) THEN
C-         Coord2(:) = Clist(:,Plist(DAT,I_P2))
C-         IF (Coord2(1) .LT. Coord1(1))THEN
C-           I_Coord = Plist(DAT,I_P1)
C-           CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
C-           Call Add_coord_after(I_P2,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-           GOTO 7900
C-         ELSE IF(Coord2(1) - Coord1(1) .LT. eps)THEN
C-           CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
C-           GOTO 7900
C-         END IF
C-           IF (debug) THEN
C-             WRITE(*,*) '8100 sort'
C-             CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
C-           END IF
C-         Coord1 = Coord2
C-         I_P1 = I_P2
C-         GOTO 8100
C-       END IF
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 8000
C-  8050 CONTINUE
C- C>
C-       IF(Debug)THEN
C-         Write(*,*) 'Scan'
C-         CALL Write_surfaces (Root_scan,Slist,Plist,Clist,RSlist)
C-         Write(*,*) 'Rootp before'
C-         CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C-       END IF
C- C
C- C ======================================================================
C- C
C- C * Start pruning by scanlines
C- C
C-       I_Scan = Root_scan
C-  9400 CONTINUE
C-       IF(I_Scan .EQ. Nullp) GOTO 9450
C-       I_PScan=Slist(DAT,I_Scan)
C-       IF (I_PScan .EQ. Nullp) THEN
C-         GOTO 9450
C-       END IF
C-  9500 CONTINUE
C-       Xmin = Clist(1,Plist(DAT,I_PScan))
C- C
C-       CALL Remove_duplicates (Rootp,
C-      &  Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
C- C
C-       I_SurfP = Rootp
C-       CALL Erase_list(Root_out,
C-      &  Slistp,Slist,Plistp,Plist,Clistp,M_Coord,M_Pnt,M_Surf)
C-  9000 CONTINUE
C-       IF(I_SurfP .EQ. Nullp) GOTO 9050
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF (I_P1 .EQ. Nullp) THEN
C-         I_SurfP = Slist(SUCC,I_SurfP)
C-         GOTO 9000
C-       END IF
C-       Coord1(:) = Clist(:,Plist(DAT,I_P1))
C-  9100 CONTINUE
C-       I_P2=Plist(SUCC,I_P1)
C-       IF (I_P2 .EQ. Nullp) THEN
C-          Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
C-       ELSE
C-          Coord2(:) = Clist(:,Plist(DAT,I_P2))
C-       ENDIF
C- C cut Left
C-       IF( Abs(Coord1(1)-Coord2(1)) .GT. Eps ) THEN
C-         Lmin = Coord1(1) .LE. Coord2(1)
C-         R=(Xmin-Coord1(1))/(coord2(1)-Coord1(1))
C- CC      write(*,*) I_P1, Coord1
C- CC      write(*,*) I_P2, Coord2
C- CC      Write(*,*) R, Xmin
C-         IF (R .GE. 0.0 .AND. R .LE. 1.0) THEN
C-           Coord1(1)=Xmin+eps
C-           Coord1(2)=R*(Coord2(2)-Coord1(2))+Coord1(2)
C-           IF(Lmin)THEN
C-              I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
C-              Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-              Coord(:)=Coord1(:)-(/2*eps,0.0D0/)
C-              I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
C-              Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-           ELSE
C-              Coord(:)=Coord1(:)-(/2*eps,0.0D0/)
C-              I_Coord = New_coord(M_Coord,Clistp,Clist,Coord)
C-              Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-              I_Coord = New_coord(M_Coord,Clistp,Clist,Coord1)
C-              Call Add_coord_after(I_P1,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-           END IF
C-           IF (debug) THEN
C-             WRITE(*,*) '9100 cut left'
C-             CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C-           END IF
C-         END IF       
C-       END IF
C-       Coord1(:)=Coord2(:)
C-       I_P1 = I_P2
C-       IF (I_P1 .NE. Nullp) GOTO 9100
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 9000
C-  9050 CONTINUE
C- C
C- C ----------------------------------------------------------
C- C
C- C * Cut points to the left of current Xmin
C- C
C-       I_SurfP = Rootp
C-  9200 CONTINUE
C-       I_Surfout = Nullp
C-       IF(I_SurfP .EQ. Nullp) GOTO 9250
C-       I_P1=Slist(DAT,I_SurfP)
C-  9300 CONTINUE
C-       IF (I_P1 .EQ. Nullp) THEN
C-         I_SurfP = Slist(SUCC,I_SurfP)
C-         GOTO 9200
C-       END IF
C-       I_Coord = Plist(DAT,I_P1)
C-       Coord1(:) = Clist(:,I_Coord)
C- 
C-       IF( Coord1(1) .LT. Xmin ) THEN
C- C        CALL Remove_coord(I_Coord,M_Coord,Clistp)
C-         I_PS = Plist(PRE,I_P1)
C-         CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
C-         IF( I_Surfout .EQ. Nullp) THEN
C-           I_Surfout = New_surf(M_Surf,Root_out,Slistp,Slist,RSlist,
C-      &        RSlist(TRANS,I_SurfP))
C-         END IF
C-         Call Add_coord(I_Surfout,I_Coord,Slist,M_Pnt,Plistp,Plist)
C-         IF (debug) THEN
C-           WRITE(*,*) '9300 removed left'
C-           CALL Write_surfaces (Root_out,Slist,Plist,Clist,RSlist)
C-         END IF
C-         IF( I_PS .NE. Nullp)THEN
C-           I_P1 = I_PS
C-         ELSE
C-           I_P1=Slist(DAT,I_SurfP)
C-           GOTO 9300
C-         END IF
C-       ELSE
C- C Keep the point
C-       END IF
C-       I_P1 = Plist(SUCC,I_P1)
C-       IF (I_P1 .NE. Nullp) GOTO 9300
C-       IF(I_Surfout .NE. Nullp) THEN
C-         CALL Remove_duplicates (I_Surfout,
C-      &    Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
C-           IF(Debug)THEN
C-             Write(*,*) 'removed dupl root_out'
C-             CALL Write_surfaces (Root_out,Slist,Plist,Clist,RSlist)
C-           END IF
C-         CALL Compute_area (I_Surfout,Slist,Plist,Clist,RSlist,eps)
C-       END IF
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 9200
C-  9250 CONTINUE
C- C>
C-       IF(Debug)THEN
C-         Write(*,*) 'root_out'
C-         CALL Write_surfaces (Root_out,Slist,Plist,Clist,RSlist)
C-         Write(*,*) 'rootp'
C-         CALL Write_surfaces (Rootp,Slist,Plist,Clist,RSlist)
C-       END IF
C- C
C- C ======================================================================
C- C
C- C * Sort by area
C- C
C-       CALL Sort_by_area (Root_out,  Slist, RSlist)
C-       Trans_w_size = Trans_w_size +
C-      &    Trans_surfaces(Root_out,Slist,Plist,Clist,RSlist,eps)
C- C
C-       IF(Debug)THEN
C-         Write(*,*) "Trans_w_size:",Trans_w_size
C-       END IF
C-       I_PScan=Plist(SUCC,I_PScan)
C-       IF (I_PScan .NE. Nullp) THEN
C-         GOTO 9500
C-       END IF
C-  9450 CONTINUE
C-       IF(Debug)THEN
C-         Write(*,*) Trans_w_size, 1.0 - Trans_w_size / Total_area
C-       END IF
C-       Transp_result = 1.0 - Trans_w_size / Total_area
C- C
C-       RETURN
C-       END

C
C ======================================================================
C ======================================================================
C
C * New_surf
C
      INTEGER FUNCTION New_surf(M_Surf,Rootp,Slistp,Slist,RSlist,Transp)
      INTEGER Nullp, SUCC, PRE, DAT, TRANS
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3, TRANS=1)
      INTEGER Slistp(*),Slist(3,*),M_Surf,Rootp,Ip
      DOUBLE PRECISION Transp, RSlist(2,*)
      M_Surf=M_Surf+1
      Ip=Slistp(M_Surf)
      IF( Rootp .NE. Nullp) THEN
        Slist(PRE,Rootp) = Ip
      ENDIF
      Slist(DAT,Ip) = NullP
      RSlist(TRANS,Ip) = Transp
      Slist(SUCC,Ip) = Rootp
      Slist(PRE,Ip) = Nullp
      Rootp = Ip
      New_surf = Ip
      RETURN
      END
C
C ======================================================================
C
C * Remove_surf
C
      SUBROUTINE Remove_surf(ISurf,Rootp,M_Surf,Slistp,Slist)
      INTEGER Nullp, SUCC, PRE
      PARAMETER (Nullp=-1, SUCC=1, PRE=2)
      INTEGER Slistp(*),Slist(3,*),M_Surf,ISurf,Rootp
      INTEGER ISs
C
      IF( Rootp .EQ. ISurf)THEN
         Rootp = Slist(SUCC,ISurf)
      END IF
      Slistp(M_Surf)=ISurf
      M_Surf=M_Surf-1
      ISs = Slist(SUCC,ISurf)
      IF (ISs .NE. Nullp) THEN
        Slist(PRE,ISs)=Slist(PRE,ISurf)
      END IF
      ISs = Slist(PRE,ISurf)
      IF (ISs .NE. Nullp) THEN
        Slist(SUCC,ISs)=Slist(SUCC,ISurf)
      END IF
      RETURN
      END
C
C ======================================================================
C
C * Add_coord
C
      SUBROUTINE Add_coord(I_Surf,I_Coord,Slist,M_Pnt,Plistp,Plist)
      INTEGER Nullp, SUCC, PRE, DAT
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER I_Surf, Slist(3,*),Plistp(*),Plist(3,*),I_Coord,M_Pnt
      INTEGER Ip, Ip1
      M_Pnt=M_Pnt+1
      Ip=Plistp(M_Pnt)
      Ip1=Slist(DAT,I_Surf)
      IF( Ip1 .NE. Nullp) THEN
        Plist(PRE,Ip1) = Ip
      ENDIF
      IF(Ip .EQ. Nullp) Write(*,*) 'Error in add_coord',M_Pnt
      Slist(DAT,I_Surf)=Ip
      Plist(PRE,Ip)=Nullp
      Plist(SUCC,Ip)=Ip1
      Plist(DAT,IP)=I_Coord
      RETURN
      END
C
C ======================================================================
C
C * Add_coord_after
C
      SUBROUTINE Add_coord_after(IC,I_Coord,Slist,M_Pnt,Plistp,Plist)
      INTEGER Nullp, SUCC, PRE, DAT
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER Slist(3,*),Plistp(*),Plist(3,*), IC, I_Coord, M_Pnt
      INTEGER Ip, Ips
      IF(IC .EQ. Nullp) THEN
        write(*,*) 'Error Null in add_coord_after'
      END IF
      M_Pnt=M_Pnt+1
      Ip=Plistp(M_Pnt)
      Ips=Plist(SUCC,IC)
      Plist(SUCC,Ip)=Ips
      Plist(PRE,Ip)=IC
      IF( Ips .NE. Nullp) THEN
         Plist(PRE,Ips)=Ip
      END IF
      Plist(SUCC,IC)=Ip
      Plist(DAT,IP)=I_Coord
      RETURN
      END
C
C ======================================================================
C
C * Remove_point
C
      SUBROUTINE Remove_point(IP,IS,M_Pnt,Plistp,Plist,Slist)
      INTEGER Nullp, SUCC, PRE, DAT
      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER Plistp(*),Plist(3,*),M_Pnt,Ip,IS,Slist(3,*)
      INTEGER Ips
      IF(IP .EQ. Nullp) THEN
        write(*,*) 'ERROR Remove_point'
        return
      END IF
      IF( Slist(DAT,IS) .EQ. IP)THEN
         Slist(DAT,IS) = Plist(SUCC,Ip)
      END IF
      Plistp(M_Pnt)=IP
      M_Pnt=M_Pnt-1
      Ips = Plist(SUCC,Ip)
      IF (Ips .NE. Nullp) THEN
        Plist(PRE,Ips)=Plist(PRE,Ip)
      END IF
      Ips = Plist(PRE,Ip)
      IF (Ips .NE. Nullp) THEN
        Plist(SUCC,Ips)=Plist(SUCC,Ip)
      END IF
      RETURN
      END
C
C ======================================================================
C
C * Remove_coord
C
      SUBROUTINE Remove_coord(IC,M_Coord,Clistp)
C      INTEGER Nullp, SUCC, PRE, DAT
C      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      INTEGER Clistp(*),M_Coord,IC
      Clistp(M_Coord)=IC
      M_Coord=M_Coord-1
      RETURN
      END
C
C ======================================================================
C
C * New_coord
C
      INTEGER FUNCTION New_coord(M_Coord,Clistp,Clist,Coords)
C      INTEGER Nullp, SUCC, PRE, DAT
C      PARAMETER (Nullp=-1, SUCC=1, PRE=2, DAT=3)
      DOUBLE PRECISION Coords(2),Clist(2,*)
      INTEGER Clistp(*),M_Coord,Ip
      M_Coord=M_Coord+1
      Ip=Clistp(M_Coord)
      Clist(:,Ip)=Coords(:)
      New_coord = Ip
      RETURN
      END

C ======================================================================
C
C * dist
C
      DOUBLE PRECISION FUNCTION dist(Coord1, Coord2)
      DOUBLE PRECISION Coord1(2), Coord2(2)
      DOUBLE PRECISION x, y
C
      X=Coord1(1)-Coord2(1)
      Y=Coord1(2)-Coord2(2)
C
      dist = sqrt( x*x + y*y )
      RETURN
      END
      

C ======================================================================
C
C * solve2
C
      LOGICAL FUNCTION  solve2 (A,X,B)
C
C Solve a 2x2 equation system A*X=B
C
C * Inputs
C    A                - Left hand side
C    B                - Right hand side
C
C * Outputs:
C     solve2        - .TRUE. if successfull
C    X                - solution
C
      DOUBLE PRECISION A(2,2), X(2), B(2)
      DOUBLE PRECISION epsilon, SLASK, R
      PARAMETER (epsilon = 1D-10)
      INTEGER I
      LOGICAL LPIV
C
      LPIV = .FALSE.
C
      IF( ABS(A(1,1)) .LE. epsilon)THEN
        IF( ABS(A(2,1)) .LE. epsilon)THEN
          solve2 = .FALSE.
          RETURN
        ELSE
          LPIV = .TRUE.
        END IF
      END IF
C
      IF( ABS(A(1,1)) .LE. ABS(A(2,1)) )THEN
        LPIV = .TRUE.
      END IF
      IF (LPIV) THEN
C Row Pivot
        DO I = 1, 2
             SLASK = A(1,I)
          A(1,I) = A(2,I)
          A(2,I) = SLASK
        END DO
        SLASK = B(1)
        B(1) = B(2)
        B(2) = SLASK
      END IF
C
      R = A(2,1) / A(1,1)
C Not used      A(2,1) = A(2,1) - A(1,1) * R (= 0.0)
      A(2,2) = A(2,2) - A(1,2) * R
      B(2) = B(2) - B(1) * R
      IF( ABS(A(2,2)) .LE. epsilon)THEN
        solve2 = .FALSE.
        RETURN
      END IF
      R = A(1,2) / A(2,2)
C Not used      A(1,2) = A(1,2) - A(2,2) * R (= 0.0)
      B(1) = B(1) - B(2) * R
      X(1) = B(1) / A(1,1)
      X(2) = B(2) / A(2,2)
      solve2 = .TRUE.
      RETURN      
      END
C
C ======================================================================
C
C * intersect
C
      LOGICAL FUNCTION intersect (P1, P2, Q1, Q2, S)
C
C Compute the intersection of the line between P1 and P2 and the
C line Q1 and Q2.
C
C * Inputs
C    P1,P2            - Points of the first line
C    Q1,Q2            - Points of the second line
C
C * Outputs:
C     intersect    - .TRUE. if successfull
C    S                - The point of intersection
C
      DOUBLE PRECISION P1(2), P2(2), Q1(2), Q2(2), S(2)
      DOUBLE PRECISION A(2,2), X(2), B(2), eps
      PARAMETER (eps = 1.0D-8)
      LOGICAL solve2, Debug
C
      Debug = .FALSE.
C
      A(1,1) = P2(1) - P1(1)
      A(1,2) = Q1(1) - Q2(1)
      A(2,1) = P2(2) - P1(2)
      A(2,2) = Q1(2) - Q2(2)
      B(1) = Q1(1) - P1(1)
      B(2) = Q1(2) - P1(2)
      IF (solve2(A,X,B)) THEN
         IF (X(1) .GE. +eps .AND. X(1) .LE. 1.0D0-eps .AND.
     &       X(2) .GE. +eps .AND. X(2) .LE. 1.0D0-eps) THEN
           intersect = .TRUE.
           S(1) = X(1) * (P2(1) - P1(1)) + P1(1)
           S(2) = X(1) * (P2(2) - P1(2)) + P1(2)
           IF(Debug)THEN
             Write(*,*) P1(1), S(1), P2(1)
             Write(*,*) P1(2), S(2), P2(2)
             Write(*,*) Q1(1), X(2) * (Q2(1) - Q1(1)) + Q1(1), Q2(1)
             Write(*,*) Q1(2), X(2) * (Q2(2) - Q1(2)) + Q1(2), Q2(2)
           END IF
         ELSE
           intersect = .FALSE.
         END IF
      ELSE
         intersect = .FALSE.
      ENDIF
      RETURN
      END
C
C ======================================================================
C
C * Remove_duplicates (Now not used [PG 20/02/2002]
C
C-       SUBROUTINE Remove_duplicates (Rootp,
C-      &  Slistp,Slist,Plistp,Plist,Clist,Clistp,M_Coord,M_Pnt,eps)
C-       INTEGER Nullp, SUCC, DAT
C-       PARAMETER (Nullp=-1, SUCC=1, DAT=3)
C-       INTEGER Slistp(*),Slist(3,*),Plistp(*),Plist(3,*),Clistp(3,*)
C-       INTEGER Rootp,M_Coord,M_Pnt
C-       DOUBLE PRECISION Clist(2,*),eps, dist
C- C
C- C Remove duplicate points in all surfaces in the list Rootp
C- C
C-       INTEGER I_SurfP, I_P1, I_P2
C-       DOUBLE PRECISION Coord1(2), Coord2(2)
C- C
C-       I_SurfP = Rootp
C-  4000 CONTINUE
C-       IF(I_SurfP .EQ. Nullp) GOTO 4050
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF (I_P1 .EQ. Nullp) THEN
C-         GOTO 4300
C-       END IF 
C-       Coord1(:) = Clist(:,Plist(DAT,I_P1))
C-  4100 CONTINUE
C-       I_P2=Plist(SUCC,I_P1)
C-       IF (I_P2 .EQ. Nullp) THEN
C-          Coord2(:) = Clist(:,Plist(DAT,Slist(DAT,I_SurfP)))
C-       ELSE
C-          Coord2(:) = Clist(:,Plist(DAT,I_P2))
C-       ENDIF
C- C cut Top and Down
C-       IF( 2*eps .GT. dist(Coord1, Coord2) ) THEN
C-         CALL Remove_coord(Plist(DAT,I_P1),M_Coord,Clistp)
C-         CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
C-       END IF
C-       Coord1(:)=Coord2(:)
C-       I_P1 = I_P2
C-       IF (I_P1 .NE. Nullp) GOTO 4100
C-  4300 CONTINUE
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 4000
C-  4050 CONTINUE
C- C
C-       RETURN      
C-       END

C
C ======================================================================
C
C * Compute_area	 (not used [PG 20/02/2002])
C
C-       SUBROUTINE Compute_area (Rootp,
C-      &  Slist,Plist,Clist,RSlist,eps)
C-       INTEGER Nullp, SUCC, DAT, AREA
C-       PARAMETER (Nullp=-1, SUCC=1, DAT=3, AREA=2)
C-       INTEGER Slist(3,*),Plist(3,*)
C-       INTEGER Rootp
C-       DOUBLE PRECISION Clist(2,*),RSlist(2,*),eps
C- C
C- C Remove duplicate points in all surfaces integer the list Rootp
C- C
C-       INTEGER I_SurfP, I_P1, I_P2, Ic, I
C- cle>  DOUBLE PRECISION Coord(2,20), Xmin, Xmax, Ymin, Ymax
C-       DOUBLE PRECISION Coord(2,100), Xmin, Xmax, Ymin, Ymax
C- C
C- cle>
C-       INTEGER SurfP
C- 
C-       DOUBLE PRECISION AB(2), AC(2), ASurf, ATriangle,
C-      &  Coord1(2), Coord2(2)
C- c     DOUBLE PRECISION ATriangle0
C- 
C-       SurfP = Rootp
C- C     Fetch 1st point A
C-       I_P1=Slist(DAT,SurfP)
C-       IF(I_P1 == NullP) THEN
C- c       CalcArea = 0
C-         RSlist(AREA,SurfP) = 0.0
C-         RETURN
C-       END IF
C-       Coord1(:) = Clist(:,Plist(DAT,I_P1))
C- 
C- C     Fetch 2nd point B
C-       I_P1=Plist(SUCC,I_P1)
C-       IF(I_P1 == NullP) THEN
C- c       CalcArea = 0
C-         RSlist(AREA,SurfP) = 0.0
C-         RETURN
C-       END IF
C-       Coord2(:) = Clist(:,Plist(DAT,I_P1))
C- c     WRITE(*,*) ' Compute_area: '
C- c     WRITE(*,*) ' Coord1: ',Coord1
C- c     WRITE(*,*) ' Coord2: ',Coord2
C-       AB = Coord2 - Coord1
C- 
C- C Repeat from here for each C point (nr 3,..)
C- c     ATriangle0 = 0.0
C-       ASurf = 0.0
C-  7200 CONTINUE
C-       I_P1=Plist(SUCC,I_P1)
C-       IF(I_P1 /= NullP) THEN
C-         Coord2(:) = Clist(:,Plist(DAT,I_P1))
C- c       WRITE(*,*) ' Coord2: ',Coord2
C-         AC = Coord2 - Coord1
C-         ATriangle = 0.5 * (AB(1)*AC(2) - AC(1)*AB(2))
C- c       if (ATriangle*ATriangle0 < 0.0) then
C- c         WRITE(*,*) ' Compute_area: Surf not convex'
C- c       end if
C- c       ATriangle0 = Atriangle
C-         ASurf = ASurf + Atriangle
C- 
C- C Next point
C-         AB = AC
C-         GOTO 7200
C-       END IF
C- 
C-       RSlist(AREA,SurfP) = ABS(ASurf)
C-       RETURN
C- cle>
C-       Xmin=1.0D10
C-       Xmax=-1.0D10
C-       Ymin=1.0D10
C-       Ymax=-1.0D10
C- C
C-       I_SurfP = Rootp
C-  4000 CONTINUE
C-       Ic=0
C-       IF(I_SurfP .EQ. Nullp) GOTO 4050
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF (I_P1 .EQ. Nullp) THEN
C-         GOTO 4300
C-       END IF 
C-       Ic=Ic+1
C-       Coord(:,Ic) = Clist(:,Plist(DAT,I_P1))
C-  4100 CONTINUE
C-       I_P2=Plist(SUCC,I_P1)
C-       IF (I_P2 .NE. Nullp) THEN
C-         Ic=Ic+1
C-         Coord(:,Ic) = Clist(:,Plist(DAT,I_P2))
C-       ENDIF
C-       I_P1 = I_P2
C-       IF (I_P1 .NE. Nullp) GOTO 4100
C-  4300 CONTINUE
C- CC        write(*,*) 'Ic = ',Ic      
C- CDeb    if (ic .gt. 20) then
C- CDeb      write(*,*) ' Compute_area: Rootp = ',Rootp
C- CDeb      WRITE(*,*) ' RSlist(AREA,I_SurfP) ',RSlist(AREA,I_SurfP)
C- CDeb      write(*,*) ' Compute_area: Ic = ',Ic
C- CDeb    end if
C-       IF (Ic .LE. 2)THEN
C-         RSlist(AREA,I_SurfP) = 0.0D0
C-       ELSE IF (Ic .EQ. 3)THEN
C-         DO I = 1, 3
C-           Xmin=Min(Xmin,Coord(1,I))
C-           Xmax=Max(Xmax,Coord(1,I))
C-           Ymin=Min(Ymin,Coord(2,I))
C-           Ymax=Max(Ymax,Coord(2,I))
C-         END DO
C- CC        write(*,*) xmin,xmax,ymin,ymax
C-         RSlist(AREA,I_SurfP) = (Xmax-Xmin)*(Ymax-Ymin)/2.0
C-       ELSE IF (Ic .EQ. 4)THEN
C-         DO I = 1, 3
C-           Xmin=Min(Xmin,Coord(1,I))
C-           Xmax=Max(Xmax,Coord(1,I))
C-           Ymin=Min(Ymin,Coord(2,I))
C-           Ymax=Max(Ymax,Coord(2,I))
C-         END DO
C- CC        write(*,*) xmin,xmax,ymin,ymax
C-         RSlist(AREA,I_SurfP) = (Xmax-Xmin)*(Ymax-Ymin)/2.0
C-         Xmin=1.0D10
C-         Xmax=-1.0D10
C-         Ymin=1.0D10
C-         Ymax=-1.0D10
C-         DO I = 3, 4
C-           Xmin=Min(Xmin,Coord(1,I))
C-           Xmax=Max(Xmax,Coord(1,I))
C-           Ymin=Min(Ymin,Coord(2,I))
C-           Ymax=Max(Ymax,Coord(2,I))
C-         END DO
C-         I=1
C-         Xmin=Min(Xmin,Coord(1,I))
C-         Xmax=Max(Xmax,Coord(1,I))
C-         Ymin=Min(Ymin,Coord(2,I))
C-         Ymax=Max(Ymax,Coord(2,I))
C- CC        write(*,*) xmin,xmax,ymin,ymax
C-         RSlist(AREA,I_SurfP) = RSlist(AREA,I_SurfP) +
C-      &     (Xmax-Xmin)*(Ymax-Ymin)/2.0
C-       ELSE
C- CC        Write(*,*) 'Error Ic',Ic
C-         DO I = 1, Ic
C- CC           Write(*,*) Coord(:,I)
C-         END DO
C-       END IF
C- C      I_SurfP = Slist(SUCC,I_SurfP)
C- C      GOTO 4000
C-  4050 CONTINUE
C- C
C- CDeb      if (ic .gt. 20) then
C- CDeb      WRITE(*,*) ' RSlist(AREA,I_SurfP) ',RSlist(AREA,I_SurfP)
C- CDeb      end if
C-       RETURN
C-       END
C
C ======================================================================
C
C * Write_surfaces (for debugging)
C
      SUBROUTINE Write_surfaces (Rootp,
     &  Slist,Plist,Clist,RSlist)
      INTEGER Nullp, SUCC, DAT, AREA
      PARAMETER (Nullp=-1, SUCC=1, DAT=3, AREA=2)
      INTEGER Slist(3,*),Plist(3,*)
      INTEGER Rootp
      DOUBLE PRECISION Clist(2,*),RSlist(2,*)
C
      INTEGER I_SurfP, I_P1
      DOUBLE PRECISION Coord1(2)
C
      I_SurfP = Rootp
 7200 CONTINUE
      IF(I_SurfP .EQ. Nullp) GOTO 7250
      Write(*,*) 'Area:',RSlist(AREA,I_SurfP)
      I_P1=Slist(DAT,I_SurfP)
      if(i_p1 .ne. nullp) then
        Coord1(:) = Clist(:,Plist(DAT,I_P1))
        Write(*,*) I_SurfP, I_P1, Coord1
      end if
 7300 CONTINUE
      IF (I_P1 .NE. Nullp) THEN
         I_P1=Plist(SUCC,I_P1)
         IF (I_P1 .NE. Nullp) THEN
            Coord1(:) = Clist(:,Plist(DAT,I_P1))
            Write(*,*) I_SurfP, I_P1, Coord1
            GOTO 7300
         END IF
      END IF
      I_SurfP = Slist(SUCC,I_SurfP)
      GOTO 7200
 7250 CONTINUE
      RETURN
      END
C
C ======================================================================
C
C * Erase_list	(Not used - [PG 20/02/2002])
C
C-       SUBROUTINE Erase_list (Rootp,
C-      &  Slistp,Slist,Plistp,Plist,Clistp,M_Coord,M_Pnt,M_Surf)
C-       INTEGER Nullp, SUCC, DAT
C-       PARAMETER (Nullp=-1, SUCC=1, DAT=3)
C-       INTEGER Slistp(*),Slist(3,*),Plistp(*),Plist(3,*),Clistp(3,*)
C-       INTEGER Rootp,M_Coord,M_Pnt,M_Surf
C- C
C-       INTEGER I_SurfP, I_P1, I_P2
C- C
C-       I_SurfP = Rootp
C-  7200 CONTINUE
C-       IF(I_SurfP .EQ. Nullp) GOTO 7250
C-       I_P1=Slist(DAT,I_SurfP)
C-  7300 CONTINUE
C-       IF (I_P1 .NE. Nullp) THEN
C-         I_P2=Plist(SUCC,I_P1)
C-         CALL Remove_coord(Plist(DAT,I_P1),M_Coord,Clistp)
C-         CALL Remove_point(I_P1,I_SurfP,M_Pnt,Plistp,Plist,Slist)
C-         I_P1=I_P2
C-         GOTO 7300
C-       END IF
C-       CALL Remove_surf(I_SurfP,Rootp,M_Surf,Slistp,Slist)
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 7200
C-  7250 CONTINUE
C-       RETURN
C-       END
C
C ======================================================================
C
C * Sort_by_area (Not used - [PG 20/02/2002])
C
C-       SUBROUTINE Sort_by_area (Rootp,  Slist, RSlist)
C-       INTEGER Nullp, SUCC, PRE, AREA
C-       PARAMETER (Nullp=-1, SUCC=1, PRE=2, AREA=2)
C-       INTEGER Rootp, Slist(3,*)
C-       DOUBLE PRECISION RSlist(2,*)
C- C
C-       INTEGER I_S1, I_S2
C-  7900 CONTINUE
C-       I_S1 = Rootp
C-  8000 CONTINUE
C-       IF(I_S1 .EQ. Nullp) GOTO 8050
C-       I_S2 = Slist(SUCC,I_S1)
C-       IF(I_S2 .EQ. Nullp) GOTO 8050
C-       IF( RSlist(AREA,I_S1) .GT. RSlist(AREA,I_S2))THEN
C-         IF( I_S1 .EQ. Rootp)THEN
C-           Rootp = I_S2
C-         END IF
C-         IF( Slist(SUCC,I_S2) .NE. Nullp) THEN
C-           Slist(PRE,Slist(SUCC,I_S2)) = I_S1
C-           Slist(SUCC,I_S1) = Slist(SUCC,I_S2)
C-         ELSE
C-           Slist(SUCC,I_S1) = Nullp
C-         END IF
C-         IF( Slist(PRE,I_S1) .NE. Nullp) THEN
C-           Slist(SUCC,Slist(PRE,I_S1)) = I_S2
C-           Slist(PRE,I_S2) = Slist(PRE,I_S1)
C-         ELSE
C-           Slist(PRE,I_S2) = Nullp
C-         END IF
C-         Slist(SUCC,I_S2) = I_S1
C-         Slist(PRE,I_S1) = I_S2
C-         GOTO 7900
C-       END IF
C-       I_S1 = I_S2
C-       GOTO 8000
C-  8050 CONTINUE
C-       RETURN
C-       END

C
C ======================================================================
C
C * Trans_surfaces  (Not used - [PG 20/02/2002])
C
C-       DOUBLE PRECISION FUNCTION Trans_surfaces (Rootp,
C-      &  Slist,Plist,Clist,RSlist,eps)
C-       INTEGER Nullp, SUCC, DAT
C-       PARAMETER (Nullp=-1, SUCC=1, DAT=3)
C-       INTEGER Slist(3,*),Plist(3,*)
C-       INTEGER Rootp
C-       DOUBLE PRECISION Clist(2,*),RSlist(2,*),eps
C- C
C-       INTEGER I_SurfP
C-       DOUBLE PRECISION Tloc, Xmin, Xmax, Y1min, Y1max, Y2min, Y2max
C- C>ML:990108
C-       DOUBLE PRECISION Xmin1, Xmax1
C- C>
C- cle>  DOUBLE PRECISION Coord(2,20)
C-       DOUBLE PRECISION Coord(2,100)
C-       INTEGER IC, I_P1, I_P2, I, Icmax
C-       INTEGER N_points
C-       LOGICAL UNDEFX,UNDEF1,UNDEF2
C- C
C-       UNDEFX=.TRUE.
C-       UNDEF1=.TRUE.
C-       UNDEF2=.TRUE.
C-       Xmin=1.0D10
C-       Xmax=-1.0D10
C-       Y1min=1.0D10
C-       Y2min=1.0D10
C-       Y1max=-1.0D10
C-       Y2max=-1.0D10
C-       I_SurfP = Rootp
C-       Tloc = 0.0D0
C-       Icmax=0
C- C      
C-  7200 CONTINUE
C-       Ic=0
C-       IF(I_SurfP .EQ. Nullp) GOTO 7250
C-       IF(N_points(I_Surfp,Slist,Plist) .LE. 2)THEN
C-         I_SurfP = Slist(SUCC,I_SurfP)
C-         GOTO 7200
C-       END IF
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF (I_P1 .EQ. Nullp) THEN
C-         GOTO 4300
C-       END IF 
C-       Ic=Ic+1
C-       Coord(:,Ic) = Clist(:,Plist(DAT,I_P1))
C-  4100 CONTINUE
C- CC      write(*,*) 'I_P1',I_P1
C-       I_P2=Plist(SUCC,I_P1)
C-       IF (I_P2 .NE. Nullp) THEN
C-         Ic=Ic+1
C-         Coord(:,Ic) = Clist(:,Plist(DAT,I_P2))
C-       ENDIF
C-       I_P1 = I_P2
C-       IF (I_P1 .NE. Nullp) GOTO 4100
C-  4300 CONTINUE
C-       Icmax= Max(Ic,Icmax)
C- C>ML:990108
C-       Xmin1=1.0D10
C-       Xmax1=-1.0D10
C-       DO I=1,IC
C-         Xmin1=Min(Xmin1,Coord(1,I))
C-         Xmax1=Max(Xmax1,Coord(1,I))
C-       END DO
C-       IF(IC .GT. 2 .AND. Xmax1-Xmin1 > 2*eps) THEN
C- C>
C-         DO I=1,IC
C-           UNDEFX=.FALSE.
C-           Xmin=Min(Xmin,Coord(1,I))
C-           Xmax=Max(Xmax,Coord(1,I))
C-         END DO
C-         DO I=1,IC
C-           IF(Abs(Coord(1,I)-Xmin) .LE. 2*eps)THEN
C-             UNDEF1=.FALSE.
C-             Y1min=Min(Y1min,Coord(2,I))
C-             Y1max=Max(Y1max,Coord(2,I))
C-           ELSE
C-             UNDEF2=.FALSE.
C-             Y2min=Min(Y2min,Coord(2,I))
C-             Y2max=Max(Y2max,Coord(2,I))
C-           END IF
C-         END DO
C- C>ML:990108
C-       END IF
C- C>
C- C>      Tloc = Tloc +
C- CA>    & RSlist(AREA,I_SurfP) * (1.0 - RSlist(TRANS,I_SurfP))
C-       I_SurfP = Slist(SUCC,I_SurfP)
C-       GOTO 7200
C-  7250 CONTINUE
C-       IF(Icmax .LE. 2) THEN
C-         UNDEFX = .TRUE.
C-         UNDEF1 = .TRUE.
C-         UNDEF2 = .TRUE.
C-       END IF
C-       IF(.NOT. (UNDEFX .OR. UNDEF1 .OR. UNDEF2))THEN
C- CC        write(*,*) 'icmax',icmax
C- CC        write(*,*) xmax,xmin
C- CC        write(*,*) y1max,y1min
C- CC        write(*,*) y2max,y2min
C-         Tloc = (Y1Max - Y1Min) * (Xmax - Xmin) / 2.0
C-         Tloc = Tloc + (Y2Max - Y2Min) * (Xmax - Xmin) / 2.0
C-       END IF
C-       Trans_surfaces = Tloc
C-       RETURN
C-       END

C ======================================================================
C
C * N_Points  (Not used - [PG 20/02/2002])
C
C-       INTEGER FUNCTION  N_Points (I_Surfp,Slist,Plist)
C-       INTEGER Nullp, SUCC, DAT
C-       PARAMETER (Nullp=-1, SUCC=1, DAT=3)
C-       INTEGER Slist(3,*),Plist(3,*)
C- C
C-       INTEGER I_SurfP
C-       INTEGER IC, I_P1
C- C
C-  7200 CONTINUE
C-       Ic=0
C-       IF(I_SurfP .EQ. Nullp) GOTO 7250
C-       I_P1=Slist(DAT,I_SurfP)
C-       IF (I_P1 .EQ. Nullp) THEN
C-         GOTO 7250
C-       END IF 
C-  4100 CONTINUE
C-       Ic=Ic+1
C-       I_P1=Plist(SUCC,I_P1)
C-       IF (I_P1 .NE. Nullp) THEN
C-         GOTO 4100
C-       ENDIF
C-  7250 CONTINUE
C-       N_Points = Ic
C-       RETURN
C-       END

C ======================================================================
C
C * Scalar_product
C
      DOUBLE PRECISION FUNCTION Scalar_product(X,Y)
      DOUBLE PRECISION X(*),Y(*)
C
C
      Scalar_product=X(1)*Y(1)+X(2)*Y(2)+X(3)*Y(3)
      RETURN
      END

C ======================================================================
C
C * Cross_product
C
      SUBROUTINE Cross_product(res,vec1,vec2)
      DOUBLE PRECISION res(*), vec1(*), vec2(*)
      res(1)=vec1(2)*vec2(3) - vec2(2)*vec1(3)
      res(2)=vec1(3)*vec2(1) - vec2(3)*vec1(1)
      res(3)=vec1(1)*vec2(2) - vec2(1)*vec1(2)
      RETURN
      END

C ======================================================================
C
C * CutZZero
C
      SUBROUTINE CutZZero(C1,C2)
      DOUBLE PRECISION C1(3), C2(3)
      DOUBLE PRECISION R
C * Cut the line C1C2 with Z=0 Return the part that is positive Z
      IF(C1(3) .LE. 0)THEN
         IF(C2(3) .LE. 0)THEN
C     * Do nothing
         ELSE 
            R=C1(3)/(C2(3)-C1(3))
            C1=C1-R*(C2-C1)
         END IF
      ELSE IF(C2(3) .LE. 0)THEN
         R=C2(3)/(C1(3)-C2(3))
         C2=C2-R*(C1-C2)
      END IF
      RETURN
      END




**********************************************************************
      SUBROUTINE ORDER(arrin, indx, N)
*
* Returns the array INDX of orders in ARRIN:
*   INDX(i) is the position of i-th sorted value (1 - smallest)
*           in the input array ARRIN
* Uses Heap sorting algorithm

      INTEGER N, indx(N)
      DOUBLE PRECISION arrin(N)

      INTEGER l, ir, indxt, i, j
	  DOUBLE PRECISION q

      do j = 1, n
         indx(j) = j
      end do

      if (n .eq. 1) return

      l = n/2 + 1
      ir= n
  10  continue

      if (l .gt. 1) then
         l = l - 1
         indxt = indx(l)
         q = arrin(indxt)
      else
         indxt = indx(ir)
         q = arrin(indxt)
         indx(ir) = indx(1)
         ir = ir - 1
         if (ir .eq. 1) then
            indx(1) = indxt
C			goto 30
            return
         end if
      end if
      i = l
      j = l + l
  20  if (j .le. ir) then
         if (j .lt. ir) then
            if (arrin(indx(j)) .lt. arrin(indx(j+1))) j = j + 1
         end if
         if (q .lt. arrin(indx(j))) then
            indx(i) = indx(j)
            i = j
            j = j + j
         else
            j = ir + 1
         end if
         go to 20
      end if
      indx(i) = indxt
      go to 10
c 30   WRITE(*,*) (indx(i),i=1,N)
c      WRITE(*,*) (arrin(i),i=1,N)
c      WRITE(*,*) 'slut'
c      stop
      END