# = NumRu::GAnalysis::BetaPlane : A class to represent a bata plane
#
# Planetary parameters are taken from the Planet module

require "numru/gphys"
require 'numru/ganalysis/planet'

module NumRu
  module GAnalysis
    class BetaPlane
      def initialize(lat0_or_latary)
        if lat0_or_latary.respond_to?(:mean)
          @lat0 = lat0_or_latary.mean.to_f
        else
          @lat0 = lat0_or_latary
        end
        @phi0 = lat0 * Math::PI / 180.0
        @f0 = 2 * Planet.omega * Math::sin(phi0)
        @beta = 2 * Planet.omega * Math::cos(phi0) / Planet::radius
        @a = Planet::radius 
      end

      attr_reader :lat0, :phi0, :f0, :beta, :a

      # Derive the x and y from the lon and lat coordinates
      #
      # ARGUMENTS
      # * gphys (GPhys) : a GPhys object (having a lon&lat as coordinates)
      #
      # RETURN VALUE
      # * [x, y] (GPhys objects)
      def get_x_y(gphys)
        lam, phi, = Planet::get_lambda_phi(gphys)
        x = lam * (@a * Math::cos(@phi0))
        x.units = @a.units
        y = ( phi - @phi0 ) * @a
        y.units = @a.units
        [x, y]
      end

      # Horizontal gradient
      #
      # ARGUMENTS
      # * gphys (GPhys) : a GPhys object (having a lon&lat as coordinates)
      # * x [GPhys or nil] : the x coordinate, which can be obtained by
      #   the get_x_y method. If nil, internally calculated by get_x_y method.
      # * y [GPhys or nil] : the y coordinate, which can be obtained by
      #   the get_x_y method. If nil, internally calculated by get_x_y method.
      #
      # RETURN VALUE
      # * [grad_x, grad_y] (GPhys objects)
      def grad_h(gphys, x=nil, y=nil)
        lond, latd = Planet::find_lon_lat_dims(gphys)
        x, y = get_x_y(gphys) if !x || !y
        bc = GPhys::Derivative::CYCLIC_OR_LINEAR
        [ gphys.cderiv(lond,bc,x), gphys.threepoint_O2nd_deriv(latd,bc,y) ]
      end

      def div_h(u, v, x=nil, y=nil)
        lond, latd = Planet::find_lon_lat_dims(u)
        x, y = get_x_y(u) if !x || !y
        bc = GPhys::Derivative::CYCLIC_OR_LINEAR
        gx = u.cderiv(lond,bc,x)
        gy = v.threepoint_O2nd_deriv(latd,bc,y)
        div = gx + gy
        div.name = "div"
        div.long_name = "div(#{u.long_name},#{v.long_name})"
        div
      end

      # Gradient in the x direction
      #
      # ARGUMENTS
      # * gphys (GPhys) : a GPhys object (having a longitude as a coordinate)
      # * x [GPhys or nil] : the x coordinate, which can be obtained by
      #   the get_x_y method. If nil, internally calculated by get_x_y method.
      #
      # RETURN VALUE
      # * grad_x (GPhys objects)
      def grad_x(gphys, x=nil)
        lond, latd = Planet::find_lon_lat_dims
        x, = get_x_y(gphys) if !x
        gphys.cderiv(lond,x)
      end

      # Gradient in the y direction
      #
      # ARGUMENTS
      # * gphys (GPhys) : a GPhys object (having a latitude as a coordinate)
      # * y [GPhys or nil] : the y coordinate, which can be obtained by
      #   the get_x_y method. If nil, internally calculated by get_x_y method.
      #
      # RETURN VALUE
      # * grad_y (GPhys objects)
      def grad_y(gphys, y=nil)
        lond, latd = Planet::find_lon_lat_dims
        x, y = get_x_y(gphys) if !y
        gphys.threepoint_O2nd_deriv(latd,y)
      end

    end
  end
end