ヤング率その他を一発で出すスクリプト

young.rb

# young.rb - ヤング率を算出するスクリプト(最小二乗法使用)
# Author: Masayoshi Wada
# 使い方わからなかったら聞いてください -> @masawada

# 以下、長さの単位はmmで

rawa = [] # aの測定値 (3点以上)
rawb = [] # bの測定値 (3点以上)
l =  # ナイフエッジ間の長さ
ucl = 0.5 # 指定
d =  # 鏡とレーザポインタ間の距離
ucd =  # 主観で決めてよい (計算値の不確かさに大きく影響)
r =  # 鏡の指示棒と測定棒間距離
ucr = 0.05 # 指定

rawms = [ # 測定値 [S/mm, m/g]
] 

# rawmsの値を逆に入れてしまったときにスワップ
#rawms.map!{|x| [x[1], x[0]]}

# DON'T EDIT FROM HERE
g = 9.7978872 # 重力加速度 (東京)
ucg = 1/(10**9) # 重力加速度の不確かさ

# functions
# 標準不確かさ計算関数
def getAverage(data)
  sum1 = 0
  data.each{|x| sum1 += x}
  return sum1/data.size.to_f
end

def getUncertainty(data)
  sum2 = 0
  avg = getAverage(data)
  data.each{|x| sum2 += (x - avg)**2}
  return Math.sqrt(sum2/(data.size-1).to_f)
end

# 最小二乗法計算関数
def linest(data)
  n = data.size
  sx = sx2 = sy = sxy = v = 0
  data.each{|p|
    sx += p[0]
    sx2 += p[0]**2
    sy += p[1]
    sxy += p[0]*p[1]
  }

  d = (n*sx2-sx**2).to_f
  d1 = sy*sx2-sx*sxy
  d2 = n*sxy-sx*sy
  a = d1/d
  b = d2/d
  wa = d/sx2
  wb = d/n
  data.each{|p| v += ((a+p[0]*b)-p[1])**2 }
  s = Math.sqrt(v.to_f/(n-2))
  return {
    :a => a,
    :b => b,
    :uca => s/Math.sqrt(wa),
    :ucb => s/Math.sqrt(wb)
  }
end

# ヤング率計算関数
def getYoung(g, ucg, l, ucl, d, ucd, a, uca, b, ucb, r, ucr, ms, ucms)
  e = (g/2)*((l**3*d)/(a**3*b*r))*(ms)
  uce = e*Math.sqrt((ucg/g)**2+(3*ucl/l)**2+(ucd/d)**2+(3*uca/a)**2+(ucb/b)**2+(ucr/r)**2+(ucms/ms)**2)
  return {
    :e => e,
    :uce => uce
  }
end

# main
# 単位変換
rawa.map!{|x| x.to_f/1000}
rawb.map!{|x| x.to_f/1000}
l = l.to_f/1000
ucl = ucl.to_f/1000
d = d.to_f/1000
ucd = ucd.to_f/1000
r = r.to_f/1000
ucr = ucr.to_f/1000
rawms.map!{|x|
  [x[0].to_f/1000, x[1].to_f/1000]
}

# 値と不確かさの取得
# aとbの値と不確かさ
a = getAverage(rawa)
uca = getUncertainty(rawa)
b = getAverage(rawb)
ucb = getUncertainty(rawb)

# ms(グラフの傾き)の値と不確かさ
retms = linest(rawms)
ms = retms[:b]
ucms = retms[:ucb]

# ヤング率の値と不確かさ
rete = getYoung(g, ucg, l, ucl, d, ucd, a, uca, b, ucb, r, ucr, ms, ucms)
e = rete[:e]
uce = rete[:uce]

# 値の表示
print "a: ", a, " m, uca: ", uca, " m\n"
print "b: ", b, " m, ucb: ", ucb, " m\n"
print "l: ", l, " m, ucl: ", ucl, " m\n"
print "d: ", d, " m, ucd: ", ucd, " m\n"
print "r: ", r, " m, ucr: ", ucr, " m\n"
print "ms: ", ms, " kg/m, ucms: ", ucms, " kg/m\n"
print "e: ", e/10**10, " *10^10 Pa, uce: ", uce/10**10, " *10^10 Pa\n"