//*********** TEIL 0 *******************
//Konstante PI merken:
PI = 3.1415926535897932384626433832795
//Matrix mit Meßdaten laden:
M=fscanfMat("messdaten.txt");

//Matrix-Dimensionen feststellen (Zeit, Iproportional, Pulse, U):
[zeilen,spalten] = size(M)

//t = t - tmin:
tmin = M(1,1);
for i=1:zeilen
    M(i,1) = M(i,1) - tmin;
end
t = M(:,1);

//*********** TEIL 1 *******************
//Aufarbeitung des Stromverlaufs
//1. Umrechnung als Strom I: Messung U an Widerstand R=98,9Ohm, Vorzeichen drehen:
imotor = M(:,2)./(-98.9);
clf();
plot(t,imotor);
 

Code 3.1.12-1: Scilab Auswerteprogramm Teil 1

//*********** TEIL 2 *******************
//Zeitpunkte merken, wo Pulsflanke -1 unterschreitet merken:
//1. durchzählen
anzahl = 0;
i=1;
gefunden = 0;
while(i<=zeilen)
    while(M(i,3)<-1.0)
        i=i+1;
        if(gefunden==0)
            anzahl = anzahl+1;
        end
        gefunden = 1;
    end
    gefunden = 0;
    i=i+1;
end
anzahl

//Zeitpunkt der Pulsflanke
t_puls = zeros(anzahl+1,1);
//Winkel der Pulsflanke
phi_puls = zeros(anzahl+1,1);

anzahl = 1;
//Bei der um 0,002s verschobenen Funktion,
//ist der Nulldurchgang um Nullpunkt:
t_puls(1,1) = 0.0;
phi_puls(1,1) = 0.0;

i=1;
gefunden = 0;

//2. merken
while(i<=zeilen)
    while(M(i,3)<-1.0)
        i=i+1;
        if(gefunden==0)
            anzahl = anzahl+1;
//Zeit um 0,002s verschieben:
            t_puls(anzahl,1) = M(i,1)-0.002;
            phi_puls(anzahl,1) = phi_puls(anzahl-1,1) +PI/6.0;
        end
        gefunden = 1;
    end
    gefunden = 0;
    i=i+1;
end

clf();
plot(t_puls,phi_puls);

//*********** TEIL 3 *******************
//Annäherung des Verlaufs durch Polynom 8. Ordnung ohne Geradenanteil:
//y(t) = c8*t^8 + c7*t^7 + c6*t^6 + c5*t^5 + c4*t^4 + c3*t^3 + c2*t^2
//Darstellung in der Form Bc+d=0 => Spalten von B: t^2,t^3,t^4, ..., d=-y
//zur Anwendung von LSQ
s2=t_puls.*t_puls;
s3=s2.*t_puls;
s4=s3.*t_puls;
s5=s4.*t_puls;
s6=s5.*t_puls;
s7=s6.*t_puls;
s8=s7.*t_puls;
B = [s2,s3,s4,s5,s6,s7,s8];
d = -phi_puls;
c = lsq(B,-d)
//Kurve t,phi mit Funktion erstellen,
//Dabei die Zeitverschiebung um 0,002s wieder rückgängig machen:
phi = zeros(zeilen,1);
for i=1:zeilen
    tt = t(i,1)-0.002;
    yy = (((((((c(7)*tt+c(6))*tt+c(5))*tt+c(4))*tt+c(3))*tt+c(2))*tt+c(1))*tt*tt);
    if(tt<0.0)
        phi(i,1)=0.0;
    else
        phi(i,1) = yy;
    end
end

//t_puls auch wieder zurückverschieben:
t_puls = t_puls+0.002.*ones(anzahl,1);

clf();
plot(t,phi,'red',t_puls,phi_puls,'blu');

//Bestimmung der momentanen Winkelgeschwindigkeit
//als Ableitung des Näherungspolynoms:
dphi = zeros(zeilen,1);
for i=1:zeilen
    tt = t(i,1)-0.002;
    yy = ((((((8.0*c(7)*tt+7.0*c(6))*tt+6.0*c(5))*tt+5.0*c(4))*tt+4.0*c(3))*tt+3.0*c(2))*tt+2.0*c(1))*tt;
    if(tt<0.0)
        dphi(i,1)=0.0;
    else
        dphi(i,1) =yy; 
    end
end

clf();
plot(t,dphi);
 

Code 3.1.12-2: Scilab Auswerteprogramm Teil 2 und 3

//*********** TEIL 4 *******************
//Der Stromverlauf soll nun auch über Funktionen angenähert werden.
//Eine gute Näherung an die Gesamtfunktion gibt die Arcus-Tangens-Funktion.
//Mit ihr wird zuert ein grober Fit des Verlaufs von Hand gesucht.
//Anschließend wird der Verlauf mit einem additiven Polynom verfeintert.
//Dieser zweite Schritt findet wieder mit Hilfe der Methode der kleinsten Quadrate statt.
//Da gerade im anfänglichen Haftreibungsanteil wenige Meßpunkte vorliegen,
//werden diese interpoliert, um Polynome höheren Grades finden zu können.


maxanzneu = zeilen*8;
grenzzeit = t(zeilen,1);
tneu = linspace(0,grenzzeit,maxanzneu);
ineu = zeros(1,maxanzneu);
ineu = interp1(t,imotor(:,1),tneu,'linear');

tneu = tneu';
ineu = ineu';

[zeil,spal] = size(tneu)

grad = 20;
B = zeros(zeil,grad);
B(:,1) = tneu;
for i=2:1:grad
   B(:,i) = B(:,i-1).*tneu; 
end

d = -ineu(1:zeil,1)+0.029*atan(7000.0*tneu);
c = lsq(B,-d)
imotor_neu = zeros(zeilen,1);
for i=1:zeilen
    tt = t(i,1);
    imotor_neu(i,1) = 0.0;
    for k= grad:-1:1
         imotor_neu(i,1) = imotor_neu(i,1)*tt;
         imotor_neu(i,1) = imotor_neu(i,1)+c(k);
    end
    imotor_neu(i,1) = imotor_neu(i,1)*tt;
    imotor_neu(i,1) = imotor_neu(i,1) +0.029*atan(7000.0*tt);
end

clf();
//plot(t(:,1),imotor(:,1),'red',t(:,1),imotor_neu(:,1),'blu');

//Ableitung für den Strom bilden:
//arctan(x)' = 1/(1+x^2)
dimotor_neu = zeros(zeilen,1);
for i=1:zeilen
    tt = t(i,1);
    dimotor_neu(i,1) = 0.0;
    for k= grad:-1:1
         dimotor_neu(i,1) = dimotor_neu(i,1)*tt;
         dimotor_neu(i,1) = dimotor_neu(i,1)+c(k)*k;
    end
    dimotor_neu(i,1) = dimotor_neu(i,1) +0.029*7000.0/(1.0+7000.0*tt*7000.0*tt);
end
plot(t(:,1),dimotor_neu(:,1),'red');
 

Code 3.1.12-3: Scilab Auswerteprogramm Teil 4: Bestimmung des Stromverlaufs

//*********** TEIL 5 *******************
//Bestimmung einzelner Modellparameter:
//Es stehen nun die zeitlichen Verläufe für folgende Größen 
//sowohl im Haftreibungsbereich, als auch bei Gleitreibung zur Verfügung:
//Strom i: imotor_neu
//Stromänderungsgeschwindigkeit di/dt: dimotor_neu
//Winkelverlauf phi: phi
//Winkelgeschwindigkeit omega: dphi
//Spannung U: konstant 5 Volt.
//
//ELEKTRISCHES MODELL
//E-a) Bestimmung aller Parameter, ohne Verwendung direkt meßbaerer Parameter 
//E-a1) Haftreibungsbereich, dphi/dt = 0rad/s
//Elektrisches Modell: u=Ri+L(di/dt) <=> Ri+L(di/dt) - u =0
//Anzahl der Stützpunkte bis t=0,002:
anz=0;
i=1;
while(t(i,1)<=0.002)
    anz = anz+1;
    i=i+1;
end
anz = anz+1;

B = [imotor_neu(1:anz,1),dimotor_neu(1:anz,1)];
d = -5.0.*ones(anz,1);
c = lsq(B,-d)
//Ergebnis: R=119Ohm, L=24mH

//E-a2) Gleitreibungsbereich, dphi/dt > 0rad/s
//Elektrisches Modell: u=Ri+L(di/dt)+kw <=> Ri+L(di/dt)+kw - u =0
B = [imotor_neu(anz:zeilen,1),dimotor_neu(anz:zeilen,1),dphi(anz:zeilen,1)];
d = -5.0.*ones(zeilen-anz+1,1);
c = lsq(B,-d)
//Ergebnis: R=115Ohm, L=912mH, k=0,0167

//E-b) Bestimmung aller Parameter, mit Verwendung direkt meßbaerer Parameter 
//E-b1) Haftreibungsbereich, dphi/dt = 0rad/s
//Vorgabe von R=24Ohm
B = [dimotor_neu(1:anz,1)];
d = -5.0.*ones(anz,1)+24.0.*imotor_neu(1:anz,1);
c = lsq(B,-d)
//Ergebnis: L=36,7mH

//E-b2) Gleitreibungsbereich, dphi/dt > 0rad/s
//Vorgabe von L=12,78mH und R=24Ohm
B = [dphi(anz:zeilen,1)];
d = -5.0.*ones(zeilen-anz+1,1)+24.0.*imotor_neu(anz:zeilen,1)+12.78*dimotor_neu(anz:zeilen,1);
c = lsq(B,-d)
//Ergebnis: k=0,0544

//MECHANISCHES MODELL
//M-b) Bestimmung aller Parameter, mit Verwendung direkt meßbaerer Parameter 
//M-b1) Haftreibungsbereich, dphi/dt = 0rad/s
//einfach den Zeitpunkt 0,002s nehmen und hier km*i bilden:
km=0.0078;
H=km*imotor_neu(anz,1)
//Ergebnis:
//H=0,0003348Nm

//M-b) Bestimmung aller Parameter, mit Verwendung direkt meßbaerer Parameter 
//M-b2) Gleitreibungsbereich, dphi/dt > 0rad/s
//Mechanisches Modell: J dphi^2/dt^2 = kmi - D dphi/dt - G
//<=> J dphi^2/dt^2 - kmi + D dphi/dt + G = 0
//plot(t,ddphi);
//Verwendung des Bereichs 0,002 bis 0,005,
//in dem die Winkelbeschleunigung monoton fällt.
//Vorgabe von J=0,00000051 kgm^2
anz2=0;
i=1;
while(t(i,1)<=0.005)
    anz2 = anz2+1;
    i=i+1;
end
anz2 = anz2+1;
//anz2-anz+1=32
J=0.00000051;
B = [-1.0.*imotor_neu(anz:anz2,1),dphi(anz:anz2,1),ones(anz2-anz+1,1)];
d = J.*ddphi(anz:anz2,1);
c = lsq(B,-d)
//Ergebnis:
//km=0,0007966
//D=0,0000207 Nms
//G=-0,00075 Nm

//Vorgabe von J=0,00000051 kgm^2 und km=0,000335
km=0.0078;
B = [dphi(anz:anz2,1),ones(anz2-anz+1,1)];
d = J.*ddphi(anz:anz2,1)-km.*imotor_neu(anz:anz2,1);
c = lsq(B,-d)
//Ergebnis:
//D=0,0000207 Nms
//G=-0,00077 Nm
 

//Bereich erweitert:
//Vorgabe von J=0,00000051 kgm^2 und km=0,000335
km=0.0078;
B = [dphi(anz2:zeilen,1),ones(zeilen-anz2+1,1)];
d = J.*ddphi(anz2:zeilen,1)-km.*imotor_neu(anz2:zeilen,1);
c = lsq(B,-d)
//Ergebnis:
//D=0,0000207 Nms
//G=-0,00077 Nm
 

Code 3.1.12-4: Scilab Auswerteprogramm Teil 5: Bestimmung der Modellparameter