clear;
clc;
clearglobal;

//Parameter festlegen
b=0.3;             //[m]
m=10.0;            //[kg]
J=(2.0/3.0)*m*b*b; //[kg*m^2]
D=0.1;             //[Nm*s/rad]

ak=[1;2];
bk=[1;2];
ck=[0;0];

function ck = kreuz(ak,bk)
    ck =  ak(1,1)*bk(2,1)-ak(2,1)*bk(1,1);
endfunction

global b,m,J,D,kreuz;

y = [0,0,0,0,0,0]';

//function f = rechteSeite(t,y)
    xx     = y(1,1);
    vx     = y(2,1);
    yy     = y(3,1);
    vy     = y(4,1);
    phi    = y(5,1);
    omega  = y(6,1);

    beta = 0.1; //beliebig gewählt hier.
    bomega = 0; //Winkelgeschw. von beta.

    FAskal = 4; //[N], später unsere Stellgröße

    //Hilfsfunktionen
    TT  = [cos(phi),-sin(phi);sin(phi),cos(phi)];
    TTR = [cos(beta),-sin(beta);sin(beta),cos(beta)];

    //Besondere Behandlung von beta für betragsm. kleine Werte:
    BETAHILF = beta;
    KLEIN    = 0.000001;
    VORZEICHEN = 1;
    if beta<0 then VORZEICHEN=-1,end
    if abs(beta)<KLEIN then BETAHILF=VORZEICHEN*KLEIN,end

    //Ausgezeichnete Punkte vorausberechnen:

    

    oS_ = [0;0];
    oH_ = [-b;0];
    oP_ = [b;0];

    oS = [xx;yy];
    oH = oS + TT*[-b;0];
    oP = oS + TT*[b;0];

    //Momentanpol oQ

    oQ_ = [-b;2*b*cos(beta)/sin(BETAHILF)];
    oQ  = oS + TT*oQ_;

    //berechenbare Kraefte bestimmen:

    //Antriebskraft:
    FA_ = TTR*[FAskal;0];    
    FA  = TT*FA_;

    //Reibkräfte
    FRH_ = -D*[2*b*omega*cos(beta)/sin(BETAHILF);0];
    FRP_ = -D*[2*b*omega*cos(beta)/sin(BETAHILF);2*b*omega];

    //...im Inertialk.:
    FRH = TT*FRH_;
    FRP = TT*FRP_;

    //Zwangsbedingungen herleiten:
    CBH = cos(beta)/sin(BETAHILF);
    OMS = omega*sin(phi);
    OMC = omega*cos(phi);
    DCBH = -bomega/(sin(BETAHILF)*sin(BETAHILF));

    //[ddot x;ddot y] = ddot phi * AAA + BBB        

    AAA = b*[-sin(phi)+2*cos(phi)*CBH;cos(phi)+2*sin(phi)*CBH];    
    BBB = b*omega*[-OMC-2*OMS*CBH + 2*cos(phi)*DCBH;-OMS+2*OMC*CBH + 2*sin(phi)*DCBH];
//alfa==ddot phi
    //Bestimmung der Zwangskräfte:
    //  [I]      [ddot x;ddot y] = FH + FP + FA + FRH + FRP
    //=>[II]     alfa * AAA + BBB = FH + FP + FA + FRH + FRP
    //
    //Euler:
    //  J*alfa=  rSH'xFH'  + rSP'x FP' + rSP'x (FA' + FRP' )
    //<=> alfa=  (rSH'xFH')/J  + (rSP'x FP')/J + (rSP'x (FA' + FRP' ))/J
    //<=> alfa=  (rSH'xFH')/J  + (rSP'x FP')/J + CCC
    //<=> [III] alfa=  ( FHskal*(-b) )/J  + ( FPskal*(b*cos(beta)) )/J + CCC
    //eingesetzt in [II] die schon ersetzte Gl. der Newton-Gl. =>
// [IV] ( ( FHskal*(-b) )/J  + ( FPskal*(b*cos(beta)) )/J + CCC ) * AAA + BBB = FH + FP + FA + FRH + FRP
// <=> FHskal*AAA*(-b/J) + FPskal*AAA*(b*cos(beta)/J) + CCC*AAA + BBB = FH + FP + FA + FRH + FRP

//<=> FHskal*AAA*(-b/J) + FPskal*AAA*(b*cos(beta)/J) -FH -FP = -CCC*AAA - BBB + FA + FRH + FRP

//FH = TT*[0;FHskal] = [cos(phi),-sin(phi);sin(phi),cos(phi)]*[0;FHskal]
//FP = TT*[-FPskal*sin(beta);FPskal*cos(beta)]

    AAAAA = [0,0;0,0];

    AAAAA(1,1) = AAA(1,1)*(-b/J) + sin(phi);
    AAAAA(1,2) = AAA(1,1)*(b*cos(beta)/J) + sin(beta)*cos(phi) +cos(beta)*sin(phi);

    AAAAA(2,1) = AAA(2,1)*(-b/J)            - cos(phi);
    AAAAA(2,2) = AAA(2,1)*(b*cos(beta)/J)   + sin(beta)*sin(phi) - cos(beta)*cos(phi);


    CCC = (  kreuz([b;0],(FA_ + FRP_))  )/J;

    bbbbb =  -CCC*AAA - BBB + FA + FRH + FRP;

    //Zangskraft bestimmen:
    Fzwang = inv(AAAAA)*bbbbb;

    FHskal = Fzwang(1,1);
    FPskal = Fzwang(2,1);

    FH_ = [0;FHskal];
    FP_ = FPskal*[-sin(beta);cos(beta)];

    FH  = TT*FH_;
    FP  = TT*FP_;

    f(1,1) = vx;
    f(2,1) = (FA(1,1)+FH(1,1)+FP(1,1)+FRH(1,1)+FRP(1,1))/m;
    f(3,1) = vy;
    f(4,1) = (FA(2,1)+FH(2,1)+FP(2,1)+FRH(2,1)+FRP(2,1))/m;
    f(5,1) = omega;
    f(6,1) = (  kreuz([-b;0],FH_) + kreuz([b;0],(FA_+FP_+FRP_))  )/J;

//endfunction

Code 0-1: Aktueller Stand zum erweiterten Simulationsmodell.