% Harmonisch erregte Schwingung mit Daempfung

m = 150 ; F0 = 0.001 ; c = 30 ; k = 0.25e2 ; Omega = 10 ; tend = 60 ; xmax = 1.5e-7 ;
m = 150 ; F0 = 0.001 ; c = 12e5 ; k = 0.25e3 ; Omega = 10 ; tend =  6 ; xmax = 2e-9 ;

t      = 0 : tend/1000 : tend ;         % Auswertung fuer 1000 Zeitschritte

omega  = sqrt(c/m) ;
eta    = Omega/omega ;
D      = k/(2*sqrt(m*c)) ;

phi    = atan((2*D*eta)/(1-eta^2)) ;
if (phi < 0) phi = phi+pi ; end         % phi muss im Bereich 0 .. pi liegen

omegaD = omega*sqrt(1-D^2) ;
A      = F0/(c*sqrt((1-eta^2)^2+4*D^2*eta^2)) ;
A1     = - A*cos(phi) ;
A2     = - A*(D*omega*cos(phi)+Omega*sin(phi))/omegaD

xhom   = exp(-D*omega*t).*(A1*cos(omegaD*t)+A2*sin(omegaD*t)) ;
xpart  = A*cos(Omega*t-phi) ;
x      = xhom + xpart ;

str    = ['x(t) für \omega = ' , num2str(omega) , ' und \Omega = ' , num2str(Omega)] ; 
subplot(2,1,1) ; plot (t , x) , grid on , title  (str) , axis([0 tend -xmax xmax])
subplot(2,1,2) ; plot (t , xpart , t , xhom , 'r') , grid on , 
 title  ('Blau: x_{part}(t) = x_{st}(t) ,  Rot: x_{hom}(t)') , axis([0 tend -xmax xmax])