#
# uupk06.rb
#
# $Id: uupk06.rb,v 1.1 2000/11/23 08:19:15 keiko Exp $
#

require "narray"
require "numru/dcl"

include NumRu
include Math


n = 100
m = 5
y1 = NArray.sfloat(n)
y2 = NArray.sfloat(n)
a = NArray.sfloat(m)

#-----------------------------------------------------------------------
dt = 1.0/(n-1)
pi = 3.14159
for j in 0..m-1
  jj = j*2+1
  a[j] = (-1.0)**(j+1) *2.0/(jj*pi)
end

for i in 0..n-1
  t = dt*i*2*pi
  y1[i] = a[0]*cos(t)
  y2[i] = 0.0
  for j in 0..m-1
    jj = j*2+1
    yy = a[j]*cos(jj*t)
    y2[i]  = y2[i] + yy
  end
end

#-----------------------------------------------------------------------
rundef = DCL::glrget('RUNDEF')

iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL::gropn iws

DCL::grfrm
DCL::sglset('LCLIP', true)

DCL::usspnt(rundef, y1)
DCL::usspnt(rundef, y2)
DCL::grswnd(0.0, 1.0, rundef, rundef)
DCL::uspfit
DCL::grstrf

DCL::usdaxs

DCL::uusarp(5999, 4999)
DCL::uvdif(rundef, y1, y2)
DCL::uulin(rundef, y1)
DCL::uulin(rundef, y2)

DCL::grfrm
DCL::sglset('LCLIP', true)

DCL::usspnt(y1, rundef)
DCL::usspnt(y2, rundef)
DCL::grswnd(rundef, rundef, 0.0, 1.0)
DCL::uspfit
DCL::grstrf

DCL::usdaxs

DCL::uhdif(y1, y2, rundef)
DCL::uulin(y1, rundef)
DCL::uulin(y2, rundef)

DCL::grcls