NumRu::GAnalysis::Planet

# File ../../lib/numru/ganalysis/planet.rb, line 105
def absvor_s(u,v)
  avor = vor_s(u,v) + @@phi.sin * (2*@@Ome)
  avor.long_name = "Absolute Vorticity"
  avor.name = "avor"
  avor
end

# File ../../lib/numru/ganalysis/planet.rb, line 66
def ave_s(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  xmean = s.mean(lond)
  cos_phi = phi.cos

  lond,latd = find_lon_lat_dims(xmean)  # find latd again
  wgt = cos_phi / cos_phi.sum
  (xmean * wgt).sum(latd)
end

# File ../../lib/numru/ganalysis/planet.rb, line 87
def div_s(u,v)
  lam, phi, lond, latd  = get_lambda_phi(u)
  cos_phi = phi.cos
  du_dlam = u.cderiv(lond,@@lonbc,lam)
  dvc_dphi  = (v*cos_phi).cderiv(latd,latbc(phi),phi)
  rot = (du_dlam + dvc_dphi) / (@@R*cos_phi)
  rot.long_name = "div(#{u.name},#{v.name})"
  rot.name = "div"
  rot
end

# File ../../lib/numru/ganalysis/planet.rb, line 221
def find_lon_lat_dims(gp, err_raise=false)
  lond = nil
  latd = nil
  (0...gp.rank).each do |dim|
    crd = gp.coord(dim)
    if /^degrees?_east$/ =~ crd.get_att("units")
      lond = dim
      break
    elsif ( ( /longitude/ =~ crd.long_name || 
              /^lon/ =~ crd.name ||
              (nm=crd.get_att("standard_name") && /longitude/ =~ nm ) &&
              (crd.units =~ Units["degrees_east"]) ) )
      lond = dim
      break
    end
  end
  (0...gp.rank).each do |dim|
    next if dim == lond
    crd = gp.coord(dim)
    if /^degrees?_north$/ =~ crd.get_att("units")
      latd = dim
      break
    elsif ( ( /latitude/ =~ crd.long_name || 
              /^lat/ =~ crd.name ||
              (nm=crd.get_att("standard_name") && /latitude/ =~ nm ) &&
              (crd.units =~ Units["degrees_north"]) ) )
      latd = dim
      break
    end
  end
  if err_raise
    raise("Longitude dimension was not found") if !lond
    raise("Latitude dimension was not found") if !latd
  end
  [lond,latd]
end

# File ../../lib/numru/ganalysis/planet.rb, line 188
def get_lambda_phi(gp, err_raise=true)
  lond, latd = find_lon_lat_dims(gp, err_raise)
  lam = lond && gp.axis(lond).to_gphys.convert_units("rad") # lon in rad
  lam.units = "" if lond                            # treat as non-dim
  phi = latd && gp.axis(latd).to_gphys.convert_units("rad") # lat in rad
  phi.units = "" if latd                            # treat as non-dim
  @@lam = lam
  @@phi = phi
  [lam, phi, lond, latd]
end

# File ../../lib/numru/ganalysis/planet.rb, line 112
def grad_s(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  xs = s.cderiv(lond,@@lonbc,lam) / (@@R*cos_phi)
  ys = s.cderiv(latd,latbc(phi),phi) / @@R
  xs.long_name = "xgrad(#{s.name})"
  xs.name = "xgrad"
  ys.long_name = "ygrad(#{s.name})"
  ys.name = "ygrad"
  [xs,ys]
end

# File ../../lib/numru/ganalysis/planet.rb, line 124
def grad_sx(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  xs = s.cderiv(lond,@@lonbc,lam) / (@@R*cos_phi)
  xs.long_name = "xgrad(#{s.name})"
  xs.name = "xgrad"
  xs
end

# File ../../lib/numru/ganalysis/planet.rb, line 133
def grad_sy(s)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  ys = s.cderiv(latd,latbc(phi),phi) / @@R
  ys.long_name = "ygrad(#{s.name})"
  ys.name = "ygrad"
  ys
end

# File ../../lib/numru/ganalysis/planet.rb, line 142
def grad_sy_cosphifact(s,cosphi_exponent)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  cos_phi_factor = cos_phi**cosphi_exponent
  ys = (s*cos_phi_factor).cderiv(latd,latbc(phi),phi)/cos_phi_factor / @@R
  ys.long_name = "ygrad(#{s.name})"
  ys.name = "ygrad"
  ys
end

# File ../../lib/numru/ganalysis/planet.rb, line 19
def init(planet)
  case planet
  when Earth
    @@R = UNumeric[6.37e6, "m"]
    @@Ome = UNumeric[2*Math::PI/8.64e4,"s-1"]
  else
    raise "Unsupported predefined planet: #{planet}."
  end
end

# File ../../lib/numru/ganalysis/planet.rb, line 48
      def latbc(phi)
        # not so good
        pi2 = Math::PI/2
        eps = 0.1
        xs = phi[0..1].val 
        xe = phi[-2..-1].val
        if (pi2-xs[0].abs) / (xs[0]-xs[1]).abs < eps and
           (pi2-xe[-1].abs) / (xe[-1]-xe[-2]).abs < eps 
          GPhys::Derivative::MIRROR_B
        else
          GPhys::Derivative::LINEAR_EXT
        end
        GPhys::Derivative::LINEAR_EXT
      end

# File ../../lib/numru/ganalysis/planet.rb, line 36
def omega; @@Ome; end

# File ../../lib/numru/ganalysis/planet.rb, line 34
def omega=(o);  @@Ome = o; end

# File ../../lib/numru/ganalysis/planet.rb, line 35
def radius; @@R; end

# File ../../lib/numru/ganalysis/planet.rb, line 33
def radius=(r); @@R = r; end

# File ../../lib/numru/ganalysis/planet.rb, line 76
def rot_s(u,v)
  lam, phi, lond, latd  = get_lambda_phi(u)
  cos_phi = phi.cos
  dv_dlam = v.cderiv(lond,@@lonbc,lam)
  duc_dphi  = (u*cos_phi).cderiv(latd,latbc(phi),phi)
  rot = (dv_dlam - duc_dphi) / (@@R*cos_phi)
  rot.long_name = "rot(#{u.name},#{v.name})"
  rot.name = "rot"
  rot
end

# File ../../lib/numru/ganalysis/planet.rb, line 98
def vor_s(u,v)
  vor = rot_s(u,v)
  vor.long_name = "Relative Vorticity"
  vor.name = "vor"
  vor
end

# File ../../lib/numru/ganalysis/planet.rb, line 159
def weight_cosphi(s, cos_exp, r_exp)
  lam, phi, lond, latd  = get_lambda_phi(s)
  cos_phi = phi.cos
  ys = s * (cos_phi**cos_exp * @@R**r_exp)
  ys
end

# File ../../lib/numru/ganalysis/planet.rb, line 166
def weight_sinphi(s, sin_exp, r_exp)
  lam, phi, lond, latd  = get_lambda_phi(s)
  sin_phi = phi.sin
  ys = s * (sin_phi**sin_exp * @@R**r_exp)
  ys
end

# File ../../lib/numru/ganalysis/planet.rb, line 152
def weight_tanphi(s, tan_exp, r_exp)
  lam, phi, lond, latd  = get_lambda_phi(s)
  tan_phi = phi.tan
  ys = s * (tan_phi**tan_exp * @@R**r_exp)
  ys
end