%  options = odeset('RelTol',1e-3,'MaxStep',dtsolv,'InitialStep',dtsolv/10);
% =====================================================================
% Testerregung eines nichtlinearen Schwingers mit Frequenzdurchstimmung
% in MATLAB 
% ------------- Beispiel: Nichtlineare Rueckstellkraft ----------------
% =====================================================================
% --- Definition der Parameter, auf die im Modell zugegriffen wird
% ----------------------------------------------------------------
% --- Modellparameter

%F_mech nichtlinear
%k_mems = 216, k2_mems = 5e07


function [fexc, xp, xp1]= wib_var_Uac(Uac)
%function [fexc, xp, xp1]= testDissjwibbeler_k_mems_var_Udc(k_mems)
%d = 0.35e-6;                    % Daempfungskonstante
%d = [0.5 1 2]*1e-6;             % Variation der Dämpfung
%k_mems = 180;                   % 1. Federkonstante
%e0 = 8.854e-12;
%R = 448e-6;
%Beta = 0.5*e0*pi*R^2;
%Beta =  2.7914e-18;

global delta d k0 m xF Frest Udc fexc d0 zahl i j

m = 3.3513e-10;                  % Masse 
k1_mems = 216;                   % 1.Federkonstante
k2_mems = 5e07;                  % 2. Federkonstante
d = 10e-6;                       % Dämpfung 
d0 = 2.2e-6;                     % abstand zwischen der Platten

f0 = sqrt(k1_mems/m)/2/pi;       % Frequenz der Schwingung

delta = d/2/m;                  % Daempfungsmass

xmax = 1e-6;                    % max. Resonatoramplitude 
n = 100;                        % Anzahl der Stuetzstellen fuer
                                % nichtlin. Kraft-Weg-Funktion
dxF = 2*xmax/n;                 % Stuetzstellenabstand
xF = [ -xmax : dxF : xmax ]';   % Stuetzstellen (x

k = ones(n+1,1)*k1_mems;        %
k1 = ones(n+1,1)*k2_mems;       %
k0 = k(n/2+1);                  %            
k2 = k1(n/2+1);                 %

Frest = zeros(n+1,1);           % Initialisierung der Restoringgkraft

%--------- Frest nichtlineare quadratische Term --------------%
for i=2:n+1                     % Integration der differentiellen
Frest = k(i).*xF - k1.*(xF).^2;
end
Frest = Frest - Frest(n/2+1);
% ------------------%% ----------------------------------------%
%--------------------- Frest lineare --------------------------%
%for i=2:n+1 % Integration der differentiellen
%Frest(i) = Frest(i-1) + k(i)*dxF;
%end 

% ------------------- Erregungsparameter-----------------------%
Udc = 0.001;
%Udc = 4;

%Uac = 0.05;
%Uac = 0.1%;
%Uac = 100;
%Uac = 0.01*f0;
%1e3*1.2675e+05

%delta = 0   1.0e+05 *0.0052 fuer fecx
%fexc = [100 105 115 120 125 133 140]*1e3;      % Erregerfrequenzen
%points=30;
points=25;
Abstand=(140-100)/points;

fexc=[100:Abstand:140]*1e3;

fexc = fexc';

% --- Erregungsfunktion als Stringvariable (notw. fuer SIM-Befehl)
aexc = ' Beta*(Udc + Uac*sin(2*pi*fexc(i)*t))^2)/(d0*(d0-x(1)))*(log(d0/x(1))) ';

%   ------------ Initialisierung der Schleife -------            %

nstep = length(fexc)                          % Anzahl der Frequenzschritte
%nparam = length(k_mems)
nparam = length(Uac)

cycles = 50;
%cycles = 200;
%--------------------------------------------------------------------------%
xic = 0; % Anfangsbedingung Weg x0
vic = 0; % Anfangsbedingung Geschwind. v0

ppp = 150; % Ausgabezeitpunkte pro Periode    50
xcol = 1; % x-Spalte des 'simout'-Arrays
vcol = 2; % v-Spalte des 'simout'-Arrays
%--------------------------------------------------------------------------%

% --- Vektoren zur Speicherung des Amplituden- und Phasenganges
xp = zeros(nstep,nparam);                      % Spitzenwerte
xp1 = zeros(nstep,nparam);                     % Spitzenwerte 1. Harmonische


% --- Schleife mit Modellaufruf
% -----------------------------
%CurrentColorArray = ['b','g','r','c','m','y','k','b','g','r','c','m','y','k','b','r','g','c','m','y','k'];


%Uac = [0.02 0.09 ]; 
Uac_ind=1:length(Uac);
zahl = zeros(size(Uac));
for j=1:nparam      % Parameter-Schleife
        
    zahl(j) = Uac(j)
     
    zahl(j)
    
% Frequenz-Schleife
    for i=1:nstep                                   
    %disp([ delta(j) fexc(i) ])                      % Fortschrittsanzeige
    disp([ delta, fexc(i) ]) 
    Tsim = cycles/fexc(i);                          % Gesamtzeit pro Frequenzschritt
    dtout = 1/(ppp*fexc(i));                        % Abstand der Ausgabezeitpunkte
    dtsolv = dtout;                                 % Max. Solver-Zeitschritt
    Tvec = [ 0 : dtout : Tsim ];                    % Ausgabezeitpunkte


% ------------ODE-Funktion --------------------%
OPTIONS = odeset('Reltol',1e-3, ...
                 'MaxStep',dtsolv, ...
                 'InitialStep',dtsolv/10 );

% Aufruf der Funktion nlinear_rest_Uac durch ode45

[t,simout] = ode45( 'nlinear_rest_Uac', Tvec, [xic vic], ...
                OPTIONS );
i
j

% --- Auswertung und Speicherung
outlen = length(t); % Anzahl der Ausgabezeitpunkte
lastcyc = simout(outlen-ppp+1:outlen,xcol); % x(t) letzte Periode
      
 
   % --- Spitzenwert
xp(i,j) = ( max(lastcyc) - min(lastcyc) ) / 2;

% --- Amplitude und Phase der 1. Harmonischen durch FFT
fftvect = fft(lastcyc)*2/ppp;
xp1(i,j) = abs(fftvect(2));
phi1(i,j) = angle( fftvect(2) )*180/pi + 90;
% --- RMS-Wert
xrms(i,j) = sqrt( sum(lastcyc.^2)/ppp );
% --- Endwerte als Anfangsbedingungen fuer naechsten Frequenzschritt
% sichern
xic = simout(outlen,xcol);
vic = simout(outlen,vcol);

   end
 %[fexc, xp, xp1]= testDissjwibbeler_k_mems_var_Uac(zahl(j)) 
end

%[fexc, xp, xp1]= testDissjwibbeler_k_mems_var_Uac(zahl(j)) 

fexc;
xp;
xp1;

figure(1)
clf
hold on 
plot(fexc, xp),'r';
plot(fexc, xp),'b';
plot(fexc, xp),'m';
grid on
xlabel('Frequenz')
ylabel('Amplitude')
%legend('Uac=0.02','Uac=0.05', 'Uac=0.09');
%title('fexc-21pts, k-216, Uac-0, cyc-50, ppp-150, Frest-wibb')
legend('Uac=0.02', 'Uac=0.09');
title('fexc-26pts, k-216, Uac-0pt001, cyc-50, ppp-150, Frest-mems, d=10e-6')
title('MEMS_erweitertes Model');

end

