disp(' ')
disp(' Milgrom_oscillator.m ')
disp(' ver 1.2  July 28, 2006 ')
disp(' by Tom Irvine  Email: tomirvine@aol.com ')
disp(' ')
disp(' This program calculates the free vibration of a  ')
disp(' single-degree-of-freedom system where ');
disp('  ')
disp(' [m/ao][d^2x/dt^2]^2 sign[d^2x/dt^2] + kx = 0')
disp('  ')
disp(' The system is subjected to an initial displacement. ')
disp(' The solution method is ode45, which is a numerical method.')
disp(' ')
%
clear omegan;
clear omegad;
clear fn;
clear damp;
clear N;
clear t;
clear x;
clear w;
clear acc;
clear v0;
clear d0;
%
disp(' ')
%v0 = input(' input initial velocity (meters/sec) ');
%disp(' ')
v0=0;
d0 = input(' input initial displacement (mm) ');
d0=d0/1000.;
%
disp(' ')
fn = input(' input natural frequency (Hz) ');
%
omegan=2*pi*fn;    %  units = (rad/sec)/sqrt(meters)
omegan2=omegan^2;
%
%   State Space Solution
%
period=1/fn;
%
dt=period/100;
sr=1/dt;
dur=20*period;
N=dur/dt;
%
% out6 = sprintf('\n dt = %f sec    dur = %f sec ',dt,dur);
% disp(out6)
%
t=linspace(0,dur,N);
t=t';
%
w(1)=d0;
w(2)=v0;
%
azero=1.2e-10; 
beta =omegan2*azero;
%
[t,w] = ode45(@Milgrom_free_function,t,w,[],[],beta);
%
dmax=max(w(:,1));
vmax=max(w(:,2));
dmin=min(w(:,1));
vmin=min(w(:,2));
%
dd=beta*(w(:,1));
%
for(i=1:N)
  acc(i)=-sign(dd(i))*sqrt( abs(dd(i)) );
end
%
amax=max(acc);
amin=min(acc);
%
out6 = sprintf('\n max disp = %11.4g mm   min disp = %11.4g mm',dmax*1000,dmin*1000);
disp(out6)
%
out6 = sprintf('\n max vel = %11.4g mm/sec   min vel = %11.4g mm/sec',vmax*1000,vmin*1000);
disp(out6)
%
out6 = sprintf('\n max acc = %11.4g m/sec^2   min acc = %11.4g m/sec^2',amax,amin);
disp(out6)
%
nzc=0;
pzc=0;
for(i=1:N-1)
%
       if(dd(i) <=0 & dd(i+1) > 0)
           pzc=pzc+1;
       end
%       
	   if(dd(i) >=0 & dd(i+1) < 0)
           nzc=nzc+1;      
       end 
end
%
out6=sprintf('\n positive zero crossings = %d ',pzc);
out7=sprintf(' negative zero crossings = %d \n',nzc);
disp(out6)
disp(out7)
%
ap = dur/pzc;
%
out6 = sprintf('\n beta = %8.4g (m/sec^4) ',beta);
disp(out6)
%
out6 = sprintf('\n Matlab period = %8.4g sec = %8.4g hr',ap,ap/3600);
disp(out6)
%
out6 = sprintf('\n Matlab fn = %8.4g Hz',1/ap);
disp(out6)
%
tt=5.98*(dmax/beta)^(1/4);
out6 = sprintf('\n Theoretical period = %8.4g sec = %8.4g hr',tt,tt/3600);
disp(out6)
%
out6 = sprintf('\n Theoretical fn = %8.4g Hz',1/tt);
disp(out6)
%
disp(' ');
disp(' plot displacement?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(1);
    plot(t,w(:,1)*1000);
    xlabel('Time (sec)'); 
    ylabel('Disp (mm)');
    grid on;    
end
%
disp(' ');
disp(' plot velocity?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(2);
    plot(t,w(:,2)*1000);
    xlabel('Time (sec)'); 
    ylabel('Vel (mm/sec)');      
    grid on;   
end
%
disp(' ');
disp(' plot acceleration?');
k=input(' 1=yes  2=no ');
%
if(k==1)
    figure(3);
    plot(t,acc );
    xlabel('Time (sec)'); 
    ylabel('Acc (m/sec^2)');      
    grid on;   
end
%
disp(' ');
disp(' plot phase portrait?');
k=input(' 1=yes  2=no ');
if(k==1)
    figure(4);
    plot(w(:,1)*1000,w(:,2)*1000);
    ylabel('Vel (mm/sec)');
    xlabel('Disp (mm)');     
    grid on;   
end
%
%  check
%
% clear q;
% for(i=1:N)
%   q(i)=acc(i)^2*sign(acc(i))+omegan2*azero*w(i,1);
% end
% plot(t,q );