module cartesian {               # helix to cartesian coordinate transformation 
contains
  subroutine h2c( inv, nn, h, ii) { 
    logical,                intent( in) :: inv   # want inverse transform?
    integer, dimension( :), intent( in) :: nn    # cartesian axes (n1,n2,n3,...)
    integer, dimension( :)              :: ii    # cartesn coords (i1,i2,i3,...)
    integer                             :: h     # shift along helix
    integer                             :: axis, n123
    n123 = 1;   if( inv) h = 1
    do axis = 1, size(  nn) { 
       if( inv)                                           # cartesian to helix
                     h = h + ( ii( axis)-1)*n123
       else                                               # helix to cartesian
                     ii( axis) = mod( ( h-1)/n123, nn( axis)) + 1
       n123 = n123 * nn( axis)
       }
    }
}