#!/usr/bin/ruby
# coding: cp932

=begin
対気モデル
Atmosphere.hの移植
=end

class InternationalStandardAtmosphere
  P0 = 101325 # Pressure at sea level [N/m^2]
  RHO0 = 1.225 # Density at sea level [kg/m^3]
  T0 = 288.15 # Temperature at sea level [K]
  A0 = 340.294 # Speed of sound at sea level [m/s]
  G0 = 9.80665 # Acceleration of gravity at sea level [m/s^2]
  R = 287.04 # Real gas constant for air [m^2 / K s^2]
  GAMMA = 1.4 # adiabatic index @see http://en.wikipedia.org/wiki/Speed_of_sound
  L = -0.0065 #

  def InternationalStandardAtmosphere.temperature(h, deltaT = 0)
    return T0 + (L * h) + deltaT
  end
  
  def InternationalStandardAtmosphere.pressure(h)
    return P0 * ((1.0 + L * h / T0) ** (-G0 / (L * R)))
  end
  
  def InternationalStandardAtmosphere.density2(h, t)
    return pressure(h) / (R * t);
  end

  def InternationalStandardAtmosphere.density(h, deltaT = 0)
    return InternationalStandardAtmosphere.density2(h, temperature(h, deltaT));
  end

  def InternationalStandardAtmosphere.mach1(h, deltaT = 0)
    return Math::sqrt(GAMMA * R * temperature(h, deltaT))
  end
  
  def InternationalStandardAtmosphere.tas(staticT, p_static, p_total)
    return Math::sqrt((2.0 * GAMMA * R * staticT / (GAMMA - 1)) \
        * ((p_total / p_static)**((GAMMA - 1) / GAMMA) - 1))
  end
  
  def InternationalStandardAtmosphere.eas(p_static, p_total)
    return InternationalStandardAtmosphere.tas(T0 / P0 * p_static, p_static, p_total)
  end
  
  def InternationalStandardAtmosphere.cas(p_impact)
    return InternationalStandardAtmosphere.tas(T0, P0, p_impact + P0)
  end
  
  def InternationalStandardAtmosphere.tas2eas(tas, staticT, p_static)
    return tas * Math::sqrt((p_static / staticT) / (P0 / T0))
  end
  
  def InternationalStandardAtmosphere.eas2cas(eas, p_static)
    p_total = p_static * ((eas**2 / ((2.0 * GAMMA * R * T0) / (GAMMA - 1) * p_static / P0) + 1)**(GAMMA / (GAMMA - 1)))
    return InternationalStandardAtmosphere.cas(p_total - p_static)
  end
  
  def InternationalStandardAtmosphere.cas2eas(cas, p_static)
    p_impact = P0 * (((cas**2 / ((2.0 * GAMMA * R * T0) / (GAMMA - 1)) + 1)**(GAMMA / (GAMMA - 1))) - 1)
    return InternationalStandardAtmosphere.eas(p_static, p_static + p_impact)
  end
  
  def InternationalStandardAtmosphere.eas2tas(eas, staticT, p_static)
    return eas * Math::sqrt((P0 / T0) / (p_static / staticT))
  end
  
  def InternationalStandardAtmosphere.p_static(tas, totalT, p_total)
    return p_total * ((1.0 - (GAMMA - 1) / GAMMA * (tas ** 2) / (R * totalT * 2)) ** (GAMMA / (GAMMA - 1)))
  end
  
  def InternationalStandardAtmosphere.p_total(u, staticT, p_static)
    return p_static * ((1.0 + (GAMMA - 1) / GAMMA * (u ** 2) / (R * staticT * 2)) ** (GAMMA / (GAMMA - 1)))
  end
  
  def InternationalStandardAtmosphere.p_impact(cas)
    return p_total(cas, T0, P0) - P0
  end
  
  def InternationalStandardAtmosphere.p_static2(tas, totalT, p_impact)
    return p_impact / ((1.0 - (GAMMA - 1) / GAMMA * (tas ** 2) / (R * totalT * 2)) ** (-GAMMA / (GAMMA - 1)) - 1.0)
  end
end


if $0 == __FILE__ then
  $: << File::dirname(__FILE__)
  require 'unit.rb'
  
  h = 0
  kt_ranges = []
  mach_ranges = []
  ARGV.reject!{|arg|
    res = true
    case arg
    when /^(\d+\.?\d*)kt$/i
      kt_ranges << $1.to_f
    when /^M(?:ach|ACH)?(\d+\.?\d*|\.\d+)$/
      mach_ranges << $1.to_f
    when /^([+-]?\d+\.?\d*)ft$/i
      h = Unit::ft2m($1.to_f)
    when /^([+-]?\d+\.?\d*)m?$/i
      h = h_str.to_f
    else
      res = false
    end
    res
  }
  raise "Some argument is incorrect: #{ARGV}" unless ARGV.empty?
  
  std_tmp = InternationalStandardAtmosphere::temperature(h)
  std_p = InternationalStandardAtmosphere::pressure(h)
  puts "Alt: #{h} [m] ( approx. #{(Unit::m2ft(h) + 0.5).to_i} [ft] )"
  puts "T(#{h}) = #{std_tmp} [K]"
  puts "p(#{h}) = #{std_p} [Pa]"
  puts "rho(#{h}) = #{InternationalStandardAtmosphere::density(h)} [kg/m^3]"
  mach1 = InternationalStandardAtmosphere::mach1(h)
  puts "M1(#{h}) = #{mach1} [m/s]"
  
  kt_ranges = [100, 200, 300, 400] if kt_ranges.empty? and mach_ranges.empty?
  
  puts [:TAS_KT, :EAS_KT, :CAS_KT, :TAS, :EAS, :CAS, :Mach].collect{|v| sprintf('%-6s', v)}.join(', ') unless kt_ranges.empty?
  kt_ranges.each{|tas_kt| # kt
    tas = Unit::kt2ms(tas_kt)
    eas = InternationalStandardAtmosphere.tas2eas(tas, std_tmp, std_p)
    eas_kt = Unit::ms2kt(eas)
    cas = InternationalStandardAtmosphere.eas2cas(eas, std_p)
    cas_kt = Unit::ms2kt(cas)
    puts [tas_kt, eas_kt, cas_kt, tas, eas, cas, tas / mach1].collect{|v| sprintf('%6.2f', v)}.join(', ')
  }
  
  puts [:CAS_KT, :EAS_KT, :TAS_KT, :CAS, :EAS, :TAS, :Mach].collect{|v| sprintf('%-6s', v)}.join(', ') unless kt_ranges.empty?
  kt_ranges.each{|cas_kt| # kt
    cas = Unit::kt2ms(cas_kt)
    eas = InternationalStandardAtmosphere.cas2eas(cas, std_p)
    eas_kt = Unit::ms2kt(eas)
    tas = InternationalStandardAtmosphere.eas2tas(eas, std_tmp, std_p)
    tas_kt = Unit::ms2kt(tas)
    puts [cas_kt, eas_kt, tas_kt, cas, eas, tas, tas / mach1].collect{|v| sprintf('%6.2f', v)}.join(', ')
  }
  
  puts [:TAS_KT, :EAS_KT, :CAS_KT, :TAS, :EAS, :CAS, :Mach].collect{|v| sprintf('%-6s', v)}.join(', ') unless mach_ranges.empty?
  mach_ranges.each{|mach|
    tas = mach1 * mach
    tas_kt = Unit::ms2kt(tas)
    eas = InternationalStandardAtmosphere.tas2eas(tas, std_tmp, std_p)
    eas_kt = Unit::ms2kt(eas)
    cas = InternationalStandardAtmosphere.eas2cas(eas, std_p)
    cas_kt = Unit::ms2kt(cas)
    puts [tas_kt, eas_kt, cas_kt, tas, eas, cas, mach].collect{|v| sprintf('%6.2f', v)}.join(', ')
  }
end