#!/usr/local/bin/ruby -d

include Math


def a(n,amp,t,y_peak)
  amp*t/(PI*n)**2*(1.0/y_peak-1.0/(y_peak-t))*(cos(2*PI*n*y_peak/t)-1)
end

def b(n,amp,t,y_peak)
  amp*t/(PI*n)**2*(1.0/y_peak-1.0/(y_peak-t))*sin(2*PI*n*y_peak/t)  
end


amp=21.0
#y_peak=2.0*PI/2.0
#t=2.0*PI/15
#y_peak=t*0.99
k_jet = 15
t=2.0*PI/k_jet
y_peak=t*0.17
beta=320

zeta = Array.new(1000)


print "#y, zeta(y)\n"

N=512 #格子点数
K=170

N.times do |i|
  y_i = 2.0*PI/N*i #i番目のy座標値
  zeta[i] = 0.0    #i番目のy座標における渦度の値
  K.times do |j|
    k=j+1

    zeta[i] += -2*PI*k/t*a(k,amp,t,y_peak)*sin(2*PI*k*y_i/t)
    zeta[i] +=  2*PI*k/t*b(k,amp,t,y_peak)*cos(2*PI*k*y_i/t)

=begin
#U(y)
    zeta[i] += -t/(2.0*PI*k)*a(k,amp,t,y_peak)*sin(2*PI*k*y_i/t)
    zeta[i] +=  t/(2.0*PI*k)*b(k,amp,t,y_peak)*cos(2*PI*k*y_i/t)
=end
=begin
#zeta(y)
    zeta[i] += a(k,amp,t,y_peak)*cos(2*PI*k*y_i/t)
    zeta[i] += b(k,amp,t,y_peak)*sin(2*PI*k*y_i/t)
=end
  end

  print "#{y_i}, #{zeta[i]} \n"
end