RungeKutta.cpp

gehe zur Dokumentation dieser Datei

/*
Autor: $Author: kunkel $ State: $State: Exp $
Datum: $Date: 2005/05/30 12:35:25 $
Version: $Revision: 1.1 $
*/

#include "RungeKutta.h"

#include "Message.h"

void RungeKutta::save(ostream & FILE) {
    FILE << "!!BEGININTEGRATOR!!" << endl;
    FILE << "0 stepsize: " <<  h << endl;
    FILE << "1000 !!END INTEGRATOR!!"<< endl;
}

void RungeKutta::load(istream & FILE) {
    int option_=0;
    string s;

    // lese solange in der Datei bis Daten des Integrators kommen
    // (integrator mark)
    while (! FILE.fail() ) {
        FILE >> s;
        if(s == "!!BEGININTEGRATOR!!")
            break;
    }
    if(s != "!!BEGININTEGRATOR!!")
        cerr << "WARNING INTEGRATOR MARK NOT FOUND!" << endl;
    while (! FILE.fail() ) {
        FILE >> option_ >> s;
        switch(option_) {
            // lade Schrittweite
        case 0:
            FILE >> h;
            break;
            // hiernach kommen keine Einstellungen mehr
        case 1000:
            return;
            break;
        default:
            cerr << "WARNING read unknown option " << option_ << " "
            << s << endl;
        }
    }
}

void RungeKutta::createAttributeWindow() {
    createBasisAttributeWindow();

    // tue nichts, wenn AttributeWindow fuer diesen Integrator bereits existiert
    if(editAttributes->additionalInformation != 0)
        return;

    // allokiere Speicher fuer additionalInformation (Schrittweite)
    GLUI * aw= editAttributes->attributeWindow;
    editAttributes->additionalInformation = (void *)(new float);

    // Spinner um Schrittweite einzustellen
    GLUI_Rollout * pO = aw->add_rollout("RungeKutta Options",true);
    aw->add_spinner_to_panel(pO,"Stepwidth", GLUI_SPINNER_DOUBLE,
                             (double*)(editAttributes->additionalInformation));

    // lade aktuelle Genauigkeit
    loadAttributes();

    // aktualisiere Spinner
    aw->sync_live();
}

void RungeKutta::loadAttributes() {
    *((double *)editAttributes->additionalInformation) = h;
}

void RungeKutta::saveAttributes() {
    h = *((double *)(editAttributes->additionalInformation));
}

void RungeKutta::deleteAdditionalInformation(void * additionalInformation) {
    delete((double *)(additionalInformation));
}

RungeKutta::RungeKutta():t(0),h(0.1),n(0),objnumber(0) {}

double RungeKutta::integrate(double time) {

    // bestimme Objektanzahl
    objnumber = objects->size();
    if(objnumber==0)
        return t;
    // allokiere Speicher fuer Massen und Radien
    m = new double[objnumber] -1;
    radius = new double[objnumber] -1;
    // berechne Anzahl der DGLn
    n = 6*objnumber;
    // allokiere Speicher fuer DGLn
    double *y = new double[n] - 1;
    double *dydx = new double[n] -1;
    double *yout = new double[n] -1;

    int k=1, i=1, j=objnumber*3+1;
    // hole alte Radien, Massen, Positionen und Geschwindigkeiten
    for (objectList::iterator g=objects->begin() ; g != objects->end();g ++) {
        Object & obj = **g;
        m[k]=obj.mass;
        radius[k]=obj.radius;
        y[i]=obj.pos.x;
        y[i+1]=obj.pos.y;
        y[i+2]=obj.pos.z;
        y[j]=obj.v.x;
        y[j+1]=obj.v.y;
        y[j+2]=obj.v.z;
        k+=1;
        i+=3;
        j+=3;
    }

    // bisher berechnete Zeitspanne
    double t_ = 0;

    // solange integrieren, bis Zeitpunkt time erreicht ist
    while(t_ < time) {

        // berechne zum ersten Mal die Ableitungen dydx
        derivs(t_,y,dydx);

        // starte Integrator, wenn Kollision stattfindet, while-Schleife
        // abbrechen; falls Schritteweite h über time hinausfuehren wuerde, nur
        // mit h=time-t_ integrieren
        if(t_+h <= time) {
            t_ += rk4(y,dydx,n,t_,h,yout);
            if(!collisions.empty())
                break;
        } else {
            t_ += rk4(y,dydx,n,t_,time-t_,yout);
            if(!collisions.empty())
                break;
        }

        // kopierere die Ergebnisse
        for (int i=1; i <= n; i++) {
            y[i] = yout[i];
        }
    }
    // aktualisiere t fuer return
    t += t_;

    i=1;
    j=objnumber*3+1;
    // kopiere neue Daten in die Objektliste
    for (objectList::iterator g=objects->begin() ; g != objects->end();g ++) {
        Object & obj = **g;
        obj.pos.x=y[i];
        obj.pos.y=y[i+1];
        obj.pos.z=y[i+2];
        obj.v.x=y[j];
        obj.v.y=y[j+1];
        obj.v.z=y[j+2];
        i+=3;
        j+=3;
    }

    // gebe integrierte Zeitspanne zurueck
    return t;
}

list<Collisionspair>* RungeKutta::getCollisionslist() {
    return &collisions;
}

double RungeKutta::distance(const int i, const int j, const double x[]) {
    if(i==j)
        return 0;
    double sum=0;
    // berechne euklidischen Abstand
    for(int d=1; d<=3; d++) {
        sum += sqr(x[(i-1)*3+d]-x[(j-1)*3+d]);
    }
    return sqrt(sum);
}

double RungeKutta::acceleration(const int i, const int komponent, const double x[] ) {
    double sum=0;
    // iteriere ueber alle anderen Objekte und bereche einzeln den gravitativen
    // Einfluss
    for(int j=1; j <= objnumber;j++) {
        if (j==i)
            continue;
        // Formel fuer Gravitationskraft
        sum+=(x[(j-1)*3+komponent]-x[(i-1)*3+komponent])*m[j]
             / cubic(distance(i,j,x));
    }
    // Gravitationskonstante Gamma nicht vergessen
    return sum*GAMMA;
}

void RungeKutta::derivs(double a, double b[], double dbda[]) {
    //X'(t) = V(t-1)
    for(int i=1; i <= n/2; i++) {
        dbda[i]= b[i+n/2];
    }
    //V'(T) = A(t-1)
    for(int i=n/2+1; i <= n; i++) {
        dbda[i] = acceleration((i-n/2+2)/3 ,(((i-n/2+2))%3)+1, b) ;
    }
}

double RungeKutta::rk4(double y[], double dydx[], int n, double x, double h, double yout[]) {
    int i;
    double xh,hh,h6,*dym,*dyt,*yt;

    dym=vec(1,n);
    dyt=vec(1,n);
    yt=vec(1,n);
    hh=h*0.5;
    h6=h/6.0;
    xh=x+hh;

    for (i=1;i<=n;i++)
        yt[i]=y[i]+hh*dydx[i];
    derivs(xh,yt,dyt);
    for (i=1;i<=n;i++)
        yt[i]=y[i]+hh*dyt[i];
    derivs(xh,yt,dym);
    for (i=1;i<=n;i++) {
        yt[i]=y[i]+h*dym[i];
        dym[i] += dyt[i];
    }
    derivs(x+h,yt,dyt);
    for (i=1;i<=n;i++)
        yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);

    // Überprüfe ob Schritteweite verringert werden muß
    for (int i=1; i <= objnumber; i++) {
        for (int j=1;  j <= objnumber; j++) {
            if(i==j)
                continue;
            Vector newposi(yout[3*(i-1)+1], yout[3*(i-1)+2], yout[3*(i-1)+3]);
            Vector oldposi(y[3*(i-1)+1], y[3*(i-1)+2], y[3*(i-1)+3]);
            Vector newposj(yout[3*(j-1)+1], yout[3*(j-1)+2], yout[3*(j-1)+3]);
            Vector oldposj(y[3*(j-1)+1], y[3*(j-1)+2], y[3*(j-1)+3]);
            double dist = (oldposi - oldposj).length();
            double potdistred = (newposi - oldposi).length() +
                                (newposj - oldposj).length();

            // Ueberpruefe, ob Objekte mit aktueller Schrittweite aneinander
            // vorbei fliegen koennten, so dass Kollision nicht unbedingt
            // entscheidbar ist
            if (dist <= potdistred) {
                free_vec(yt,1,n);
                free_vec(dyt,1,n);
                free_vec(dym,1,n);
                Message::msg(Message::MSG,
                             "Collision possible: reducing stepsize");
                // halbiere Schrittweite, wenn sie nicht schon zu klein ist
                if(h > 1e-10)
                    return rk4(y, dydx, n, x, h/2, yout);
            }
        }
    }

    // ueberpruefe ob Ueberlappungen von Objekten (Kollisionsen) stattgefunden
    // haben
    for (int i=1; i <= objnumber; i++) {
        for (int j=i;  j <= objnumber; j++) {
            if(i==j)
                continue;
            Vector newposi(yout[3*(i-1)+1], yout[3*(i-1)+2], yout[3*(i-1)+3]);
            Vector newposj(yout[3*(j-1)+1], yout[3*(j-1)+2], yout[3*(j-1)+3]);
            if ((newposi - newposj).length() <= radius[i]+radius[j]) {
                // falls Ueberlappungen stattgefunden haben, speichere diese in
                // Kollisionsliste
                collisions.push_back(Collisionspair(i,j));
            }
        }
    }

    free_vec(yt,1,n);
    free_vec(dyt,1,n);
    free_vec(dym,1,n);

    // gebe durchgefuehrte Schrittweite zurueck
    return h;
}

void RungeKutta::free_vec(double *v, long nl, long nh) {
    free((FREE_ARG) (v+nl-NR_END));
}

double *RungeKutta::vec(long nl, long nh) {
    double *v;

    v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
    if (!v)
        throw Error("allocation failure in vec()");
    return v-nl+NR_END;
}

---