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


function [fexc, xp, xp1]= testDissjwibbeler_k_mems_var_Uac(Uac)
%function [fexc, xp, xp1]= testDissjwibbeler_k_mems_var_Udc(k_mems)


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

m = 3.3513e-10;                 % Masse 
k_mems = 216;
%k_mems = 180;                   % 1. Federkonstante
%k_mems = [195.0658 203.8140 212.5478 216 221.2895 230.0308 238.7725];
%k2_mems = 5e07;                % 2. Federkonstante
d = 10e-6;
%d = 0.35e-6;                    % D�mpfungskonstante
%d = [0.5 1 2]*1e-6;             % Variation der D�mpfung
 

d0 = 2.2e-6;


%e0 = 8.854e-12;
%R = 448e-6;
%Beta = 0.5*e0*pi*R^2;
Beta =  2.7914e-18;
f0 = sqrt(k_mems/m)/2/pi;       % Frequenz der Schwingung

delta = d/2/m;                  % D�mpfungsmass

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)*k_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 f�r nichtlineare quadratische Term --------------%
for i=2:n+1 % Integration der differentiellen
%Frest(i) = Frest(i-1) + k(i)*dxF - k1(i).*dxF.*dxF;
Frest(i) = Frest(i-1) + k(i)*dxF;
end 

Frest = Frest - Frest(n/2+1);
% ------------------%% ----------------------------------------%

% --------------- Erregungsparameter -----------------%
Udc = 5;
%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

%delta =  1.0e+10*0.0000
%fexc = [100 105 115 120 125 133 140]*126750000;

%delta = 0   1.0e+05 *0.0052 fuer fecx1
%fexc1 = [0.8 0.85 0.9:0.01:1 1.15 1.2]*1.2675e+05;

%delta = 0   1.0e+05 *0.0052 fuer fecx2
%fexc2 = [0.85 0.9:0.001:1 1.01:0.001:1.1 1.11]*1.2675e+05;

%fexc = [0.8 0.85 0.9:0.01:1 1.01:0.01:1.1 1.15 1.2].*f0;

%%delta = 0   1.0e+05*0.0052 fuer fecx3
%fexc3 = [100 102 103:0.090:105 106:0.090:110 112 115]*126750000;
fexc = fexc';

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


%              ----- 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.001 0.009 0.03 ]; 
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 nonlinear_restoring_var_k_mems durch ode45

[t,simout] = ode45( 'nonlinear_restoring_k_mems_var_Uac', Tvec, [xic vic], ...
                OPTIONS );
i
j
    end
 
end
[fexc, xp, xp1]= testDissjwibbeler_k_mems_var_Uac(zahl(j)) 
 
 
figure(1)   
hold on
%CurrentColor = CurrentColorArray(zahl(j));
%CurrentColor = CurrentColorArray(k_ind);
%plot_string=strcat('-',CurrentColor);

%plot(fexc, xp, CurrentColor)
%plot(fexc, xp1, CurrentColor)
plot(fexc, xp)
xlabel('Frequenz')
ylabel('Amplitude')

end




% --- 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);






