source: TOOLS/MOZAIC/src/MOZAIC/math.f90 @ 3918

! -*- Mode: f90 -*-
MODULE math
   !
   !> Some useful mathematical functions
   !
   USE declare
   IMPLICIT NONE
   !
   !> 2D vector type, quadruple precision
   TYPE v2d
      REAL (kind=rd) :: x !< X-coordinate
      REAL (kind=rd) :: y !< Y-coordinates
   END TYPE v2d
   !
   !> Multiply 2D vectors
   INTERFACE OPERATOR ( * )  
      MODULE PROCEDURE vect_prod
   END INTERFACE
   !
CONTAINS
   !
   ELEMENTAL FUNCTION vect_prod (a, b) 
      !> Vector product
      USE declare
      IMPLICIT NONE
      !
      REAL (kind=rl) :: vect_prod
      TYPE(v2d), INTENT (in) :: a !< First vector to multiply 
      TYPE(v2d), INTENT (in) :: b !< 2nd vector to multiply
      !
      vect_prod = a%x * b%y - a%y * b%x
      !
      RETURN
   END FUNCTION vect_prod
   !
   FUNCTION poly_surf (xa, ya )
      !> Surface of a polygone
      USE declare
      IMPLICIT NONE
      !
      REAL (kind=rl) :: poly_surf
      REAL (kind=rl), INTENT (in), DIMENSION ( :) :: xa, ya
      INTEGER (kind=il) :: jn
      REAL (kind=rd) :: zs
      !
      zs = 0.0_rd
      DO jn = 1, SIZE (xa) - 2
         zs = zs + v2d ( xa(jn+1)-xa(jn), ya(jn+1)-ya(jn)) * v2d ( xa(jn+2)-xa(jn), ya(jn+2)-ya(jn))
      END DO
      poly_surf = 0.5_rl * zs    
      !
      RETURN
   END FUNCTION poly_surf
   !
   ELEMENTAL FUNCTION dist (x1, y1, x2, y2)
      !> Distance between two points (cartesian)
      USE declare
      IMPLICIT NONE
      REAL (kind=rl) :: dist
      REAL (kind=rl), INTENT ( in) :: x1, y1, x2, y2
      dist = SQRT ( (x2-x1)**2 + (y2-y1)**2 )
      RETURN
   END FUNCTION dist
   ELEMENTAL FUNCTION ddist (x1, y1, x2, y2)
      !> Distance between two points (cartesian), quadruple precision
      USE declare
      IMPLICIT NONE
      REAL (kind=rd) :: ddist
      REAL (kind=rd), INTENT ( in) :: x1, y1, x2, y2
      ddist = SQRT ( (x2-x1)**2 + (y2-y1)**2 )
      RETURN
   END FUNCTION ddist
   !
   ELEMENTAL FUNCTION dgeodist (x1, y1, x2, y2)
      !> Distance between two points (on sphere), quadruple precision
      USE declare
      IMPLICIT NONE
      REAL (kind=rd) :: dgeodist
      REAL (kind=rd), INTENT (in) :: x1, y1, x2, y2
      REAL (kind=rd) :: zs, zrad
      zrad = ACOS ( -1.0_rd) / 180.0_rd
      zs =     SIN (zrad*REAL(y1,rd)) * SIN (zrad*REAL(y2,rd)) &
         &  +  COS (zrad*REAL(y1,rd)) * COS (zrad*REAL(y2,rd)) * COS(zrad*REAL (x2-x1,rd))
      zs = MAX (-1.0_rd, MIN (1.0_rd, zs))
      dgeodist =  REAL ( ACOS (zs), KIND=rl)
      RETURN
   END FUNCTION dgeodist
   !
!-$$   ELEMENTAL FUNCTION geodist (x1, y1, x2, y2)
!-$$      !> Distance between two points (on sphere)
!-$$      USE declare
!-$$      IMPLICIT NONE
!-$$      REAL (kind=rl) :: geodist
!-$$      REAL (kind=rl), INTENT (in) :: x1, y1, x2, y2
!-$$      geodist = REAL (dgeodist (REAL (x1,rd), REAL (y1,rd), REAL (x2,rd), REAL(y2,rd)),rl)
!-$$   END FUNCTION geodist

   ELEMENTAL FUNCTION geodist (zlon1, zlat1, zlon2, zlat2)
      !> Distance between two points (on sphere)
      USE declare
      IMPLICIT NONE
      REAL (kind=rl) :: geodist
      REAL (kind=rl), INTENT (in) :: zlon1, zlat1, zlon2, zlat2
      REAL (kind=rl) :: zs, zrad
      zrad = ACOS ( -1.0_rl) / 180.0_rl
      zs = SIN (zrad*zlat1) * SIN (zrad*zlat2) +  COS (zrad*zlat1) * COS (zrad*zlat2) * COS(zrad*(zlon2-zlon1))
      zs = MAX (-1.0_rl, MIN (1.0_rl, zs))
      geodist =  ACOS (zs)
      RETURN
   END FUNCTION geodist

!-$$   ELEMENTAL FUNCTION geodist_min (zlon1, zlat1, zlon2, zlat2, zdmax)
!-$$      !> Distance between two points (on sphere). Simplified for too far points
!-$$      USE declare
!-$$      IMPLICIT NONE
!-$$      REAL (kind=rl) :: geodist
!-$$      REAL (kind=rl), INTENT (in) :: zlon1, zlat1, zlon2, zlat2, zdmax
!-$$      REAL (kind=rl) :: zs, zrad
!-$$
!-$$      IF ( ABS (zlat2-zlat1)  .GT. zdmax * zrad ) THEN
!-$$         geodist = zdmax
!-$$         RETURN
!-$$      END IF
!-$$      
!-$$      
!-$$      zrad = ACOS ( -1.0_rl) / 180.0_rl
!-$$      zs = SIN (zrad*zlat1) * SIN (zrad*zlat2) +  COS (zrad*zlat1) * COS (zrad*zlat2) * COS(zrad*(zlon2-zlon1))
!-$$      zs = MAX (-1.0_rl, MIN (1.0_rl, zs))
!-$$      geodist =  ACOS (zs)
!-$$      RETURN
!-$$   END FUNCTION geodist_min
   
   !
   FUNCTION cartdist (zlon1, zlat1, zlon2, zlat2)
      !> Distance between two points (on sphere)
      ! Simplified computation
      USE declare
      IMPLICIT NONE
      REAL (kind=rl) :: cartdist
      REAL (kind=rl), INTENT (in) :: zlon1, zlat1, zlon2, zlat2
      REAL (kind=rl) :: x1, y1, z1, x2, y2, z2, c1, c2, dx, dy, dz, dxyz
      REAL (kind=rl) :: zrad
      zrad = ACOS ( -1.0_rl) / 180.0_rl
      c1   = COS(zrad*zlat1)
      x1   = c1*COS(zrad*zlon1)
      y1   = c1*SIN(zrad*zlon1)
      z1   = SIN(zrad*zlat1)
      c2   = COS(zrad*zlat2) 
      x2   = c2*COS(zrad*zlon2)
      y2   = c2*SIN(zrad*zlon2)
      z2   = SIN(zrad*zlat2)
      
      dx = x1-x1 ; dy = y1-y2 ; dz = z1-z2
      dxyz = dx*dx + dy*dy + dz*dz
      cartdist = SQRT ( dxyz )

   END FUNCTION cartdist
   !
   ELEMENTAL FUNCTION an_som (x0, y0, x1, y1, x2, y2)
      !> Angle au sommet d'un triangle spherique
      USE declare
      IMPLICIT NONE
      REAL (kind=rd) :: an_som
      REAL (kind=rd), INTENT ( in) :: x0, y0, x1, y1, x2, y2
      REAL (kind=rd) :: za, zb, zc, zs
      za = dgeodist (REAL(x1,rd), REAL(y1,rd), REAL(x2,rd), REAL(y2,rd))
      zb = dgeodist (REAL(x0,rd), REAL(y0,rd), REAL(x1,rd), REAL(y1,rd))
      zc = dgeodist (REAL(x0,rd), REAL(y0,rd), REAL(x2,rd), REAL(y2,rd))
      IF ( za == 0.0_rd .OR. zb == 0.0_rd .OR. zc == 0.0_rd ) THEN 
         an_som = 0.0_rd
      ELSE
         zs =  ( COS(za)-COS(zb)*COS(zc) ) / ( SIN(zb)*SIN(zc) )
         an_som = REAL (ACOS (zs))
      ENDIF
      RETURN
   END FUNCTION an_som
   !
   ELEMENTAL FUNCTION surf_sphe (x1, y1, x2, y2, x3, y3)  
      !> Surface d'un triangle spherique
      USE declare
      IMPLICIT NONE
      REAL (kind=rd) :: surf_sphe
      REAL (kind=rd), INTENT ( in) :: x1, y1, x2, y2, x3, y3
      REAL (kind=rd) :: za, zb, zc
      za = an_som (REAL(x1,rd), REAL(y1,rd), REAL(x3,rd), REAL(y3,rd), REAL(x2,rd), REAL(y2,rd))
      zb = an_som (REAL(x2,rd), REAL(y2,rd), REAL(x1,rd), REAL(y1,rd), REAL(x3,rd), REAL(y3,rd))
      zc = an_som (REAL(x3,rd), REAL(y3,rd), REAL(x2,rd), REAL(y2,rd), REAL(x1,rd), REAL(y1,rd))
      IF ( za == 0.0_rd .OR. zb == 0.0_rd .OR. zc == 0.0_rd ) THEN 
         surf_sphe = 0.0_rl
      ELSE
         surf_sphe = za + zb + zc - ACOS ( -1.0_rd)
      ENDIF
      RETURN
   END FUNCTION surf_sphe
   !
   FUNCTION surf_box (xed, yed)
      !> Surface d'une boite de modele
      USE declare
      IMPLICIT NONE
      REAL (kind=rl) :: surf_box
      REAL (kind=rl), DIMENSION(9), INTENT(in) :: xed, yed
      INTEGER (kind=il) :: je, jne
      REAL (kind=rd) :: zsurf

      jne = SIZE(xed)
      zsurf = 0.0_rl

      DO je = 1, jne-1
         zsurf = zsurf + surf_sphe (REAL(xed(jne ),rd), REAL(yed(jne ),rd), &
            &                       REAL(xed(je  ),rd), REAL(yed(je  ),rd), &
            &                       REAL(xed(je+1),rd), REAL(yed(je+1),rd) )
      END DO

      surf_box = REAL (zsurf, kind = rl)

      RETURN
   END FUNCTION surf_box
   !
END MODULE math