#
# uspk09.rb
#
# $Id: uspk09.rb,v 1.1 2000/11/22 14:36:03 keiko Exp $
#

require "narray"
require "numru/dcl"

include NumRu
include Math


n = 200
m = 5
x = NArray.sfloat(n)
y = NArray.sfloat(n)
yc = NArray.sfloat(n, m)
a = NArray.sfloat(m)

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

x = x.indgen * dt
for i in 0..n-1
  t = dt*i*2*PI
  y[i] = 0.0
  for j in 0..m-1
    jj = j*2+1
    yc[i,j] = a[j]*cos(jj*t)
    y[i]  = y[i] + yc[i,j]
  end
end

#--------------------------- 1ST PAGE ----------------------------------
iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
rundef = DCL::glrget('RUNDEF')
DCL::gropn iws

DCL::grfrm
DCL::grsvpt(0.2, 0.8, 0.2, 0.6)
DCL::usspnt(rundef, yc)

DCL::ussttl('X-AXIS', ' ', 'COMPONENTS', ' ')
DCL::usgrph(x, yc[true,0])

for j in 1..m-1
  ip = j.divmod(4)[1] + 1
  DCL::uulinz(x, yc[true, j], 1, ip)
end

#     --- NEW FIG ---
DCL::grfig
DCL::grsvpt(0.2, 0.8, 0.62, 0.82)
DCL::uzlset('LABELXB', false)

DCL::ussttl('X-AXIS', ' ', 'TOTAL', ' ')
DCL::usgrph(x, y)

#--------------------------- 2ND PAGE ----------------------------------
DCL::grfrm
DCL::grsvpt(0.2, 0.8, 0.2, 0.6)

DCL::uzfact(0.5)
DCL::uzlset('LABELXB', true)

DCL::usspnt(rundef, yc)
DCL::ussttl('X-AXIS', ' ', 'COMPONENTS', ' ')
DCL::usgrph(x, yc[true,0])

for j in 1..m-1
  ip = j.divmod(4)[1] + 1
  DCL::uulinz(x, yc[true, j], 1, ip)
end

#     --- NEW FIG ---
DCL::grfig
DCL::grsvpt(0.2, 0.8, 0.62, 0.82)
DCL::uzlset('LABELXB', false)

DCL::ussttl('X-AXIS', ' ', 'TOTAL', ' ')
DCL::usgrph(x, y)

DCL::grcls