proj.F90

!    Copyright 2011 Davide Cesari <dcesari69 at gmail dot com>
!
!    This file is part of FortranGIS.
!
!    FortranGIS is free software: you can redistribute it and/or modify
!    it under the terms of the GNU Lesser General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    FortranGIS is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Lesser General Public License for more details.
!
!    You should have received a copy of the GNU Lesser General Public
!    License along with FortranGIS.  If not, see
!    <http://www.gnu.org/licenses/>.

MODULE proj
use,INTRINSIC :: iso_c_binding
IMPLICIT NONE

TYPE,bind(c) :: pj_object
  PRIVATE
  TYPE(c_ptr) :: ptr = c_null_ptr
END TYPE pj_object

TYPE(pj_object),PARAMETER :: pj_object_null=pj_object(c_null_ptr)

TYPE,bind(c) :: pjuv_object
  REAL(kind=c_double) :: u=huge(1.0_c_double)
  REAL(kind=c_double) :: v=huge(1.0_c_double)
END TYPE pjuv_object

REAL(kind=c_double),PARAMETER :: pj_deg_to_rad=.0174532925199432958d0 
REAL(kind=c_double),PARAMETER :: pj_rad_to_deg=57.29577951308232d0 

INTERFACE pj_associated
  MODULE PROCEDURE pj_associated_object
END INTERFACE pj_associated

INTERFACE
  FUNCTION pj_init_plus(name) BIND(C,name='pj_init_plus')
  import
  CHARACTER(kind=c_char) :: name(*)
  TYPE(pj_object) :: pj_init_plus
  END FUNCTION pj_init_plus
END INTERFACE

INTERFACE
  FUNCTION pj_transform(src, dst, point_count, point_offset, x, y, z) &
   bind(c,name='pj_transform')
  import
  TYPE(pj_object),VALUE :: src
  TYPE(pj_object),VALUE :: dst
  INTEGER(kind=c_long),VALUE :: point_count
  INTEGER(kind=c_int),VALUE :: point_offset
  REAL(kind=c_double) :: x(*)
  REAL(kind=c_double) :: y(*)
  REAL(kind=c_double) :: z(*)
  INTEGER(kind=c_int) :: pj_transform
  END FUNCTION pj_transform
END INTERFACE

INTERFACE
  FUNCTION pj_datum_transform(src, dst, point_count, point_offset, x, y, z) &
   bind(c,name='pj_datum_transform')
  import
  TYPE(pj_object),VALUE :: src
  TYPE(pj_object),VALUE :: dst
  INTEGER(kind=c_long),VALUE :: point_count
  INTEGER(kind=c_int),VALUE :: point_offset
  REAL(kind=c_double) :: x(*)
  REAL(kind=c_double) :: y(*)
  REAL(kind=c_double) :: z(*)
  INTEGER(kind=c_int) :: pj_datum_transform
  END FUNCTION pj_datum_transform
END INTERFACE

INTERFACE
  FUNCTION pj_fwd(val, proj) BIND(C,name='pj_fwd')
  import
  TYPE(pjuv_object),VALUE :: val
  TYPE(pj_object),VALUE :: proj
  TYPE(pjuv_object) :: pj_fwd
  END FUNCTION pj_fwd
END INTERFACE

INTERFACE
  FUNCTION pj_inv(val, proj) BIND(C,name='pj_inv')
  import
  TYPE(pjuv_object),VALUE :: val
  TYPE(pj_object),VALUE :: proj
  TYPE(pjuv_object) :: pj_inv
  END FUNCTION pj_inv
END INTERFACE

INTERFACE
  FUNCTION pj_geocentric_to_geodetic(a, es, point_count, point_offset, x, y, z) &
   bind(c,name='pj_geocentric_to_geodetic')
  import
  REAL(kind=c_double),VALUE :: a
  REAL(kind=c_double),VALUE :: es
  INTEGER(kind=c_long),VALUE :: point_count
  INTEGER(kind=c_int),VALUE :: point_offset
  REAL(kind=c_double) :: x(*)
  REAL(kind=c_double) :: y(*)
  REAL(kind=c_double) :: z(*)
  INTEGER(kind=c_int) :: pj_geocentric_to_geodetic
  END FUNCTION pj_geocentric_to_geodetic
END INTERFACE

INTERFACE
  FUNCTION pj_geodetic_to_geocentric(a, es, point_count, point_offset, x, y, z) &
   bind(c,name='pj_geodetic_to_geocentric')
  import
  REAL(kind=c_double),VALUE :: a
  REAL(kind=c_double),VALUE :: es
  INTEGER(kind=c_long),VALUE :: point_count
  INTEGER(kind=c_int),VALUE :: point_offset
  REAL(kind=c_double) :: x(*)
  REAL(kind=c_double) :: y(*)
  REAL(kind=c_double) :: z(*)
  INTEGER(kind=c_int) :: pj_geodetic_to_geocentric
  END FUNCTION pj_geodetic_to_geocentric
END INTERFACE

INTERFACE
  FUNCTION pj_compare_datums(srcdefn, dstdefn) BIND(C,name='pj_compare_datums')
  import
  TYPE(pj_object),VALUE :: srcdefn
  TYPE(pj_object),VALUE :: dstdefn
  INTEGER(kind=c_int) :: pj_compare_datums
  END FUNCTION pj_compare_datums
END INTERFACE

INTERFACE
  FUNCTION pj_is_latlong(proj) BIND(C,name='pj_is_latlong')
  import
  TYPE(pj_object),VALUE :: proj
  INTEGER(kind=c_int) :: pj_is_latlong
  END FUNCTION pj_is_latlong
END INTERFACE

INTERFACE
  FUNCTION pj_is_geocent(proj) BIND(C,name='pj_is_geocent')
  import
  TYPE(pj_object),VALUE :: proj
  INTEGER(kind=c_int) :: pj_is_geocent
  END FUNCTION pj_is_geocent
END INTERFACE

INTERFACE
  FUNCTION pj_get_def(proj, options) BIND(C,name='pj_get_def')
  import
  TYPE(pj_object),VALUE :: proj
  INTEGER(kind=c_int) :: options
  TYPE(c_ptr) :: pj_get_def
  END FUNCTION pj_get_def
END INTERFACE
  

INTERFACE
  FUNCTION pj_latlong_from_proj(proj) BIND(C,name='pj_latlong_from_proj')
  import
  TYPE(pj_object),VALUE :: proj
  TYPE(pj_object) :: pj_latlong_from_proj
  END FUNCTION pj_latlong_from_proj
END INTERFACE

!INTERFACE
!  FUNCTION pj_apply_gridshift(c, i, point_count, point_offset, x, y, z) &
!   BIND(C,name='pj_apply_gridshift')
!  IMPORT
!  CHARACTER(kind=c_char) :: c(*)
!  INTEGER(kind=c_int),VALUE :: i
!  INTEGER(kind=c_long),VALUE :: point_count
!  INTEGER(kind=c_int),VALUE :: point_offset
!  REAL(kind=c_double) :: x(*)
!  REAL(kind=c_double) :: y(*)
!  REAL(kind=c_double) :: z(*)
!  INTEGER(kind=c_int) :: pj_apply_gridshift
!  END FUNCTION pj_apply_gridshift
!END INTERFACE
!
!INTERFACE
!  SUBROUTINE pj_deallocate_grids() BIND(C,name='pj_deallocate_grids')
!  END SUBROUTINE pj_deallocate_grids
!END INTERFACE

INTERFACE
  SUBROUTINE pj_free(proj) BIND(C,name='pj_free')
  import
  TYPE(pj_object),VALUE :: proj
  END SUBROUTINE pj_free
END INTERFACE

!void pj_pr_list(projPJ);
!void pj_set_finder( const char *(*)(const char *) );
!void pj_set_searchpath ( int count, const char **path );
!projPJ pj_init(int, char **);
!char *pj_get_def(projPJ, int);
!void *pj_malloc(size_t);
!void pj_dalloc(void *);
!char *pj_strerrno(int);
!int *pj_get_errno_ref(void);
!const char *pj_get_release(void);

CONTAINS

FUNCTION pj_associated_object(arg1, arg2) RESULT(associated_)
TYPE(pj_object),INTENT(in) :: arg1
TYPE(pj_object),INTENT(in),OPTIONAL :: arg2
LOGICAL :: associated_
IF(PRESENT(arg2)) THEN
  associated_ = c_associated(arg1%ptr, arg2%ptr)
ELSE
  associated_ = c_associated(arg1%ptr)
ENDIF
END FUNCTION pj_associated_object

! Fortran specific version of some functions

FUNCTION pj_transform_f(src, dst, x, y, z)
TYPE(pj_object),VALUE :: src
TYPE(pj_object),VALUE :: dst
REAL(kind=c_double) :: x(:)
REAL(kind=c_double) :: y(:)
REAL(kind=c_double),OPTIONAL :: z(:)
INTEGER(kind=c_int) :: pj_transform_f

REAL(kind=c_double),POINTER :: dummyz(:)

IF (PRESENT(z)) THEN
  pj_transform_f = pj_transform(src, dst, &
   int(min(SIZE(x),SIZE(y),SIZE(z)), kind=c_long), 1_c_int, x, y, z)
ELSE
  NULLIFY(dummyz)
  pj_transform_f = pj_transform(src, dst, &
   int(min(SIZE(x),SIZE(y)), kind=c_long), 1_c_int, x, y, dummyz)
ENDIF

END FUNCTION pj_transform_f


FUNCTION pj_datum_transform_f(src, dst, x, y, z)
TYPE(pj_object),VALUE :: src
TYPE(pj_object),VALUE :: dst
REAL(kind=c_double) :: x(:)
REAL(kind=c_double) :: y(:)
REAL(kind=c_double),OPTIONAL :: z(:)
INTEGER(kind=c_int) :: pj_datum_transform_f

REAL(kind=c_double),POINTER :: dummyz(:)

IF (PRESENT(z)) THEN
  pj_datum_transform_f = pj_datum_transform(src, dst, &
   int(min(SIZE(x),SIZE(y),SIZE(z)), kind=c_long), 1_c_int, x, y, z)
ELSE
  NULLIFY(dummyz)
  pj_datum_transform_f = pj_datum_transform(src, dst, &
   int(min(SIZE(x),SIZE(y)), kind=c_long), 1_c_int, x, y, dummyz)
ENDIF

END FUNCTION pj_datum_transform_f


FUNCTION pj_geocentric_to_geodetic_f(a, es, x, y, z)
REAL(kind=c_double),VALUE :: a
REAL(kind=c_double),VALUE :: es
REAL(kind=c_double) :: x(:)
REAL(kind=c_double) :: y(:)
REAL(kind=c_double) :: z(:)
INTEGER(kind=c_int) :: pj_geocentric_to_geodetic_f

pj_geocentric_to_geodetic_f = &
 pj_geocentric_to_geodetic(a, es, &
 int(min(SIZE(x),SIZE(y),SIZE(x)), kind=c_long), 1_c_int, x, y, z)

END FUNCTION pj_geocentric_to_geodetic_f


FUNCTION pj_geodetic_to_geocentric_f(a, es, x, y, z)
REAL(kind=c_double),VALUE :: a
REAL(kind=c_double),VALUE :: es
REAL(kind=c_double) :: x(:)
REAL(kind=c_double) :: y(:)
REAL(kind=c_double) :: z(:)
INTEGER(kind=c_int) :: pj_geodetic_to_geocentric_f

pj_geodetic_to_geocentric_f = &
 pj_geodetic_to_geocentric(a, es, &
 int(min(SIZE(x),SIZE(y),SIZE(x)), kind=c_long), 1_c_int, x, y, z)

END FUNCTION pj_geodetic_to_geocentric_f


END MODULE proj