Module: NumRu::GAnalysis::SphHarmonicIsoG

Defined in:
lib/numru/ganalysis/sph_harmonic_iso_g.rb

Overview

Spherical harmonics library for equally spaced lon-lat grids using SHTLIB in DCL

REQUIREMENT

  • This module can handle grid data that satisfy the following

    • dimension 0 (first dim) is longitude, dimension 1 is latitude.

    • must be equally spaced and cover the globe such as lon = 0, 2,.., 358 [,360], or -180, -170,..,170 [,180]. etc (in increasing order; cyclic extension is applied internally) and lat = -90,-88,..,90. (or 90, 98,..,-90; pole to pole)

  • Data missig is not allowed

Constant Summary collapse

@@work = 
@@im = @@jm = @@mm = nil

Class Method Summary collapse

Class Method Details

.__factor_n_m(&factor_n_m) ⇒ Object 



195
196
197
198
199
200
201
202
203
204
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 195

def __factor_n_m( &factor_n_m )
  ms = DCLExt.sht_get_m(@@mm)
  ns = DCLExt.sht_get_n(@@mm)
  len = ms.length
  f = NArray.float(len)
  for i in 0...len
    f[i] = factor_n_m.call(ns[i],ms[i])
  end
  f
end

.check_and_init(gphys, mm) ⇒ Object

Check the valdity of the lon-lat grid.

This method raises an exception if the data is not acceptable

ARGUMENTS

  • gp [GPhys] : data to check its grid

  • mm [Integer] : truncation wavenumber

RETURN VALUE

  • gphys [GPhys] : same as input, but for the 1st dim is cyclically extended if needed.

  • np2sp [true or false] : true if data is N.pole to S.pole, false if data is S.pole to N.pole

Raises:

  • (ArgumentError) 


36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 36

def check_and_init(gphys, mm)

  raise(ArgumentError, "Invalid rank (#{gphys.rank})") if gphys.rank < 2

  lon = gphys.coord(0).val
  if lon[0] > lon[-1]
    raise ArgumentError, "Longitude must be in the increasing order."
  end
  gphys = gphys.cyclic_ext(0)

  lat = gphys.coord(1).convert_units("degrees").val
  eps = 0.1
  if (lat[0]-90.0).abs < eps &&  (lat[-1]+90.0).abs < eps
    # N.pole to S.pole
    np2sp = true
  elsif (lat[0]+90.0).abs > eps ||  (lat[-1]-90.0).abs > eps
    raise ArgumentError, "Not pole to pole: #{lat[0]..lat[-1]}."
  else
    np2sp = false
  end

  nx, ny, = gphys.shape
  im = (nx-1)/2
  jm = (ny-1)/2

  if (mm+1)/2 > jm || mm+1 > im
    raise(ArgumentError,"mm=#{mm} is too big for im=#{im} & jm=#{jm}")
  end
  if @@work.nil? || @@mm != mm || @@jm != jm || @@im != im
    @@work = DCL.shtint(mm,jm,im)
    @@mm = mm
    @@jm = jm
    @@im = im
  end
  [gphys, np2sp]
end

.div_h(u, v, mm, radius = 1.0, &factor_n_m) ⇒ Object

Horizontal divergence on the sphere

ARGUMENTS

  • u,v [GPhys] : the x and y components to take divergence

  • mm [Integer] : truncation wavenumber

  • radius (optional; defaut=1.0) [Numeric(if non-dim) or UNumeric] radius of the sphere



173
174
175
176
177
178
179
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 173

def div_h(u,v, mm, radius=1.0, &factor_n_m )
  gp = xdiv(u, mm, &factor_n_m) + ydiv(v, mm, &factor_n_m)
  gp *= radius**(-1) if radius != 1.0
  gp.long_name = "div_h(#{u.name},#{v.name})"
  gp.name = "div_h"
  gp
end

.filt(gp, mm, deriv = nil, lap = nil, &factor_n_m) ⇒ Object

Spherical harmonics filter

ARGUMENTS

  • gp [GPhys] : grid data to filter (lon,lat,…: rank >= 2)

  • mm [Integer] : truncation wavenumber

  • deriv (optional) [nil(default) or Symbol] (let lambda be longitude and phi be latitude) if :xgrad (or :xdiv), applies del / cosphi dellambda, if :ygrad, applies del / delphi, if :ydiv, applies del cosphi / cosphi delphi.

  • lap (optional) [nil(default) or Integer] If 1, laplacian is taken; if -1 inverserser laplacian is taken – note that you can also explitly take laplacian by using the factor_n_m block (see below).

  • factor_n_m (optional block) spectral filter in the form of {|n,m| } that returns the factor if n and m are integers. Example {|n,m| -n*(n+1)} to take laplacian.



92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 92

def filt( gp, mm, deriv=nil, lap=nil, &factor_n_m )
  iswg2s = isws2g = 0
  case deriv
  when :xgrad, :xdiv
    isws2g = -1
  #when :xderiv
  ##(horinout, developpers memo: just multiply with cos\phi afterward)
  #  iswg2s = -1   # don't know if isws2g is better
  #  raise("Under development: need cos_phi factor")
  when :ygrad
    isws2g = 1
  when :ydiv
    iswg2s = 1
  when nil
    # do nothing
  else
    raise ArgumentError,"Unsupported operation #{deriv}"
  end
  nlon = gp.shape[0]
  gp, np2sp = check_and_init(gp, mm)   # sets @@work, @@mm, @@im & @@jm
  na = gp.val
  na = na.to_na if na.respond_to?(:to_na)   # for NArrayMiss
  shape = na.shape
  naf = NArray.float( *shape )  # output variable (filtered NArray)
  f = __factor_n_m( &factor_n_m ) if factor_n_m
  loop_for_3rd_or_greater_dim(shape){|sel|
    w,s = DCL.shtg2s(@@mm,iswg2s,na[*sel],@@work)
    s = DCL.shtlap(@@mm,lap,s) if lap
    s = s * f if factor_n_m
    s = -s if np2sp && (isws2g==1 || iswg2s==1 )  # y-reversed --> *(-1)
    w,g =  DCL.shts2g(@@mm,@@jm,@@im,isws2g,s,@@work)
    naf[*sel] = g
  }
  vaf = VArray.new(naf, gp.data, gp.name)
  gp = GPhys.new( gp.grid, vaf)
  if gp.shape[0] != nlon
    # cyclically extended --> trim it to the original shape
    gp = gp[0...nlon,false]
  end
  gp
end

.lapla_h(gp, mm, radius = 1.0, order = 1) ⇒ Object

Horizontal Laplacian on the sphere

  • gp [GPhys] : grid data (lon,lat,…: rank >= 2)

  • mm [Integer] : truncation wavenumber

  • radius (optional; defaut=1.0) [Numeric(if non-dim) or UNumeric] radius of the sphere



153
154
155
156
157
158
159
160
161
162
163
164
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 153

def lapla_h( gp, mm, radius=1.0, order=1 )
  if order==-1 || order==1
    gp = filt( gp, mm, nil, order )
  elsif order > 0
    gp = filt( gp, mm ){|n,m| (-n*(n+1))**order }
  else
    # negative --> avoid zero division
    gp = filt( gp, mm ){|n,m| n==0 ? 0 : (-n*(n+1))**order }
  end
  gp *= radius**(-2) if radius != 1.0
  gp
end

.loop_for_3rd_or_greater_dim(shape, &block) ⇒ Object

Raises:

  • (ArgumentError) 


206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 206

def loop_for_3rd_or_greater_dim(shape,&block)
  raise(ArgumentError, "block not given") if !block
  sh3 = shape[2..-1]
  rank3 = sh3.length
  csh3 = [1]
  (1...rank3).each{|d| csh3[d] = sh3[d-1]*csh3[d-1]}
  len = 1
  sh3.each{|n| len *= n}
  for i in 0...len
    sel = [true,true]
    (0...rank3).each do |d|
      sel.push( (i/csh3[d]) % sh3[d] )
    end
    block.call(sel)
  end
end

.rot_h(u, v, mm, radius = 1.0, &factor_n_m) ⇒ Object

Horizontal rotation on the sphere

  • u,v [GPhys] : the x and y components to take rotation

  • mm [Integer] : truncation wavenumber

  • radius (optional; defaut=1.0) [Numeric(if non-dim) or UNumeric] radius of the sphere



187
188
189
190
191
192
193
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 187

def rot_h(u,v, mm, radius=1.0, &factor_n_m )
  gp = xdiv(v, mm, &factor_n_m) - ydiv(u, mm, &factor_n_m)
  gp *= radius**(-1) if radius != 1.0
  gp.long_name = "rot_h(#{u.name},#{v.name})"
  gp.name = "rot_h"
  gp
end

.xgrad(gp, mm, &factor_n_m) ⇒ Object Also known as: xdiv 



134
135
136
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 134

def xgrad( gp, mm, &factor_n_m )
  filt( gp, mm, :xgrad, &factor_n_m )
end

.ydiv(gp, mm, &factor_n_m) ⇒ Object 



143
144
145
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 143

def ydiv( gp, mm, &factor_n_m )
  filt( gp, mm, :ydiv, &factor_n_m )
end

.ygrad(gp, mm, &factor_n_m) ⇒ Object 



140
141
142
 
# File 'lib/numru/ganalysis/sph_harmonic_iso_g.rb', line 140

def ygrad( gp, mm, &factor_n_m )
  filt( gp, mm, :ygrad, &factor_n_m )
end