require "narray"
require "numru/dcl"

############################################################

=begin
=module NumRu::Cartesian in cartesian.rb

==Index
* ((<module NumRu::Cartesian>))
  * ((<init>))
    * 初期化ルーチン. ((<NumRu::Cartesian>))のその他のメソッドを呼ぶ前に必ず呼ばねばならない.
  * ((<setwn>))
    * データのグリッドおよび切断波数から南北, 東西それぞれの分割数の 1/2 を求める. 
  * ((<xderiv>))
    * 経度方向の微分
  * ((<yderiv>))
    * 緯度方向の微分
  * ((<rot>))
    * 回転
  * ((<div>))
    * 発散
  * ((<xdiv>))
    * 経度成分の発散
  * ((<ydiv>))
    * 緯度成分の発散
  * ((<grad>))
    * 勾配
  * ((<xgrad>))
    * 経度方向の勾配
  * ((<ygrad>))
    * 緯度方向の勾配
  * ((<lapla>))
    * ラプラシア
  * ((<ilapla>))
    * ラプラシアンの逆関数
  * ((<jacobian>))
    * ヤコビアン
  * ((<ydiv_m0>))
    * 東西波数 0 のデータの緯度方向の発散を求める. 
  * ((<ygrad_m0>))
    * 東西波数 0 のデータの緯度方向の勾配を求める.
  * ((<yderiv_m0>))
    * 東西波数 0 のデータの緯度方向の微分を求める.

=module NumRu::Cartesian

Module functions for Operater of spherical corrdinate for NArray using RubyDCL.

球面座標上の主に微分演算に関する関数を集めたモジュール. 

データの格子点は少なくとも

  緯度方向 : 南緯 90 度, 赤道, 北緯 90 度
  経度方向 : 西経 180 度, 0, 東経 180 度(=西経 180 度)

の 3 つの格子点を持たねばならない. また格子間隔は全て等しくなければならない.


=end
############################################################

module NumRu

  module Cartesian
    include NumRu

    #<<< module data >>>
    @@mm = nil    # <切断波数(M)>
    @@jm = nil    # <南北分割数の1/2(J)>
    @@im = nil    # <東西分割数の1/2(I)>
    @@work = nil  # <shtlibの他のルーチンで用いられる作業領域>

#    module_function # delete commenout after test-finished 

    def init(mm, jm, im)
      
      # <<@@mm, @@jm, @@im が全て値を持っていて, かつ対応する引数と全て一致する場合のみ shtint を呼び出さない.>>
      if (!@@mm || @@mm != mm) ||
	  (!@@jm || @@jm != jm)  ||
	  (!@@im || @@im != im)  then
	@@mm = mm; @@jm = jm; @@im = im
	@@work = DCL.shtint(mm, jm, im)
      end
      # << @@mm, @@jm, @@im value is returned for test only.>>
      return @@work, @@mm, @@jm, @@im
    end
    
    def setwn(z, total_wn=nil)
=begin 
#--------------------------------------------------------------------------------------------
      # <<データの経度方向の両端が同じ値で無い場合, 糊代をつくる>>
      if (z[0, false] - z[-1, false]).abs.max > 0 
	ln = z.shape[0]-1
	z = z[[0..ln, 0],false]
      end
      if (z[true, 0, false].ne z[true, -1, false])
	# under construction <<データの南端, 北端が極と一致していない場合に調節するよう処理>>
      end
#--------------------------------------------------------------------------------------------
=end

      if z.shape[0]%2 == 0
	raise ArgumentError,"z.shape[0] must be odd."
      else
	im = (z.shape[0]-1)/2
      end
      if z.shape[1]%2 == 0
	raise ArgumentError,"z.shape[1] must be odd."
      else
	jm = (z.shape[1]-1)/2
      end
      if (!total_wn)
	mm = [im, jm].max
      else
	mm = total_wn  # i will be check value of total_wn; (mm =< 2jm-1 && mm =< im-1) 
      end
      return z, mm, jm, im
    end

    def xderiv(z, total_wn = nil)
      z, mm, jm, im = setwn(z, total_wn)
      * = init(mm, jm, im)
      wz, sz = DCL.shtg2s(@@mm,@@jm,@@im,0,z,@@work)   # グリッド   => スペクトル
      wz, gz = DCL.shts2g(@@mm,@@jm,@@im,1,sz,@@work)  # スペクトル => グリッド(経度微分)
      return gz
    end

    def yderiv(z, total_wn = nil)
      z, mm, jm, im = setwn(z, total_wn)
      * = init(mm, jm, im) 
      wz, sz = DCL.shtg2s(@@mm,@@jm,@@im,0,z,@@work)   # グリッド   => スペクトル
      wz, gz = DCL.shts2g(@@mm,@@jm,@@im,-1,sz,@@work) # スペクトル => グリッド(緯度微分)
      return gz
    end

    def div(zx, zy, total_wn = nil)
      #
    end

    def xdiv(zx, total_wn = nil)
      #
    end

    def ydiv(zy, total_wn = nil)
      #
    end

    def grad 
      #
    end

    def xgrad 
      #
    end

    def ygrad 
      #
    end

    def rot 
      #
    end


    def lapla 
      #
    end

    def ilapla 
      #
    end

    def jacobian 
      #
    end

    def ydiv_m0 
      #
    end

    def ygrad_m0 
      #
    end

    def yderiv_m0 
      #
    end

  end

end

############################################################

# dummyclass for unit test

class DummyCartesian
  include NumRu::Cartesian
end