function[timhist] = makesine(filname,upper,lower,freq,nhcyc)
%
% timhist=makesine(filename, upperlevel, lowerlevel, frequency, # 1/2 cycles)
%
% th=makesine('c:\bts\TestA.bin', 0.1, -0.3, 1.0, 3)
%
% Makes sine wave of arbitrary upper/lower levels, frequency & # 1/2 cycles
%
% timhist = time history result file with 1/2 cycle of zeros on front/back
% filname = result file name in '' chrs ( use full path, eg. 'C:\BTS\TestA.bin' )
% upper   = upper amplitude of sine wave
% lower   = lower amplitude of sine wave
% freq    = frequency (Hz) of sine wave
% nhcyc   = # of 1/2 cycles in waveform (1-?) NOT counting 1st and last 1/2
%           cycles
%
%

dt = 1.0/1024.0;             % sample rate is 1024Hz

halfcyc = 1.0/(2.0*freq);     % 1/2 cycle time period (sec)

ptshalf = halfcyc/dt;         % data points/half cycle

for i=1:ptshalf,             % leading zeros of 1/2 cycle
    timhist(i,1)=0.0;
end;
index=ptshalf;
omega=2.0*pi*freq;
upratio=upper/(upper-lower);
nfirst=upratio*ptshalf;
omega1st=omega/upratio;
time=0.0;
for j=1:nfirst,
    timhist(index+j,1)=upper*0.5*(1.0-cos(omega1st*time));      % create 1st half cycle from 0 to upper
    time=time+dt;
end;
index=index+j;
time=0.0;
for i=1:ptshalf,
    halfcyc(i)=0.5*(1.0-cos(omega*time));   % create unity 1/2 cycle
    time=time+dt;
end;
direction = -1;
for i=1:nhcyc,
    if direction < 0.0
        mult=lower-upper;       % amplitude of 1/2 cycle from lower to upper
        offset=upper;
    else
        offset=lower;
        mult=upper-lower;       % amplitude of 1/2 cycle from upper to lower
    end;
    for j=1:ptshalf,
        timhist(index+j,1)=mult*halfcyc(j)+offset;    % create next 1/2 cycle
    end;
    index=index+ptshalf;
    direction=-direction;
end;
if direction < 0.0
    mult=-upper;
    offset=upper;
else
    mult=-lower;
    offset=lower;
end;
lowratio=abs(mult)/(upper-lower);
nlast=lowratio*ptshalf;
omegalast=omega/lowratio;
time=0.0;
for j=1:nlast,
    timhist(index+j,1)=mult*0.5*(1.0-cos(omegalast*time))+offset;  % create last 1/2 cycle to zero
    time=time+dt;
end;
index=index+j;
for j=1:ptshalf,
    timhist(index+j,1)=0.0;    % last 1/2 cycles of zeros
end;
index=index+ptshalf;

savebin(filname, timhist, dt);  % save the binary file
hold off;
time=[0.0:dt:(index-1)*dt];
plot(time,timhist);
grid on;
title('makesine.m time history');
xlabel('Time (secs)');
ylabel('Amplitude');