RungeKuttaAdaptive.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 "RungeKuttaAdaptive.h"

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

void RungeKuttaAdaptive::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 Genauigkeitseinstellung
        case 0:
            FILE >> eps;
            break;
            // hiernach kommen keine Einstellungen mehr
        case 1000:
            return;
            break;
        default:
            cerr << "WARNING read unknown option " << option_ << " "
            << s << endl;
        }
    }
}

void RungeKuttaAdaptive::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 double);

    GLUI_Rollout * pO = aw->add_rollout("RungeKutta Adaptive Options",true);

    // Spinner um Genauigkeit (eps) einzustellen
    aw->add_spinner_to_panel(pO,"eps", GLUI_SPINNER_DOUBLE,
                             (double*)(editAttributes->additionalInformation));

    // lade aktuelle Genauigkeit
    loadAttributes();

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

void RungeKuttaAdaptive::loadAttributes() {
    *((double *)editAttributes->additionalInformation) =  eps;
}

void RungeKuttaAdaptive::saveAttributes() {
    eps = *((double *)(editAttributes->additionalInformation));
}

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

RungeKuttaAdaptive::RungeKuttaAdaptive() {
    n=0;
    t = 0;
    objnumber = 0;

    // empfohlene Genauigkeit
    eps=1e-8;
}

double RungeKuttaAdaptive::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;

    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;
    }

    // max. Schrittlaenge
    double stepsize = time;
    nextsteplength = 0;

    // 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
        if(rkqs(y,dydx,n,& t_, stepsize,eps,y,& actualsteplength,
                & nextsteplength))
            break;

        // max Schrittlaenge verringern
        stepsize = time-t_;
    }

    // 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>* RungeKuttaAdaptive::getCollisionslist() {
    return &collisions;
}

double RungeKuttaAdaptive::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 RungeKuttaAdaptive::acceleration(const int i, const int komponente, 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+komponente]-x[(i-1)*3+komponente])*m[j] /
             cubic(distance(i,j,x));
    }
    // Gravitationskonstante Gamma nicht vergessen
    return sum*GAMMA;
}

void RungeKuttaAdaptive::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);
    }
}


void RungeKuttaAdaptive::rkck(double y[], double dydx[], int n, double x, double h, double yout[], double yerr[]) {
    int i;
    static double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,
                                          b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2,
                                                                                  b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0,
                                                                                                                  b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0,
                                                                                                                                                         b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0,
                                                                                                                                                                                  c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0,
                                                                                                                                                                                                                   dc5 = -277.00/14336.0;

    double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0,
                                     dc4=c4-13525.0/55296.0,dc6=c6-0.25;
    double *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp;

    ak2=vec(1,n);
    ak3=vec(1,n);
    ak4=vec(1,n);
    ak5=vec(1,n);
    ak6=vec(1,n);
    ytemp=vec(1,n);
    for (i=1;i<=n;i++)
        ytemp[i]=y[i]+b21*h*dydx[i];
    derivs(x+a2*h,ytemp,ak2);
    for (i=1;i<=n;i++)
        ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
    derivs(x+a3*h,ytemp,ak3);
    for (i=1;i<=n;i++)
        ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
    derivs(x+a4*h,ytemp,ak4);
    for (i=1;i<=n;i++)
        ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);

    derivs(x+a5*h,ytemp,ak5);
    for (i=1;i<=n;i++)
        ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]
                         +b65*ak5[i]);
    derivs(x+a6*h,ytemp,ak6);
    for (i=1;i<=n;i++)
        yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
    for (i=1;i<=n;i++)
        yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
    free_vec(ytemp,1,n);
    free_vec(ak6,1,n);
    free_vec(ak5,1,n);
    free_vec(ak4,1,n);
    free_vec(ak3,1,n);
    free_vec(ak2,1,n);
}

bool RungeKuttaAdaptive::rkqs(double y[], double dydx[], int n, double *x, double htry, double eps, double yscal[], double *hdid, double *hnext) {

    int i;
    double errmax,h,htemp,xnew,*yerr,*ytemp;

    yerr=vec(1,n);
    ytemp=vec(1,n);
    h=htry;

    int counter=0;
    for (;;) {
        rkck(y,dydx,n,*x,h,ytemp,yerr);
        errmax=0.0;
        for (i=1;i<=n;i++)
            errmax=max(errmax,abs(yerr[i]/yscal[i]));
        errmax /= eps;
        if (errmax <= 1.0)
            break;
        htemp=SAFETY*h*(double) (pow(errmax,PSHRNK));
        h=(h >= 0.0 ? max(htemp,0.1*h) : min(htemp,0.1*h));
        xnew=(*x)+h;

        // falls der Integrator nicht von der Stelle kommt -> abbrechen
        if(xnew == *x)
            throw Error("stepsize underflow in rungekuttaadaptive (type 1)");

        // falls Integrator zu lange braucht -> abbrechen
        if(counter > 1000)
            throw Error("stepsize underflow in rungekuttaadaptive (type 2)");
        counter++;
    }
    if (errmax > ERRCON)
        *hnext=SAFETY*h*(double) pow(errmax,PGROW);
    else
        *hnext=5.0*h;
    *x += (*hdid=h);


    // 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(ytemp[3*(i-1)+1], ytemp[3*(i-1)+2], ytemp[3*(i-1)+3]);
            Vector newposj(ytemp[3*(j-1)+1], ytemp[3*(j-1)+2], ytemp[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));
            }
        }
    }

    for (i=1;i<=n;i++)
        y[i]=ytemp[i];
    free_vec(ytemp,1,n);
    free_vec(yerr,1,n);

    // geben false zurueck, wenn keine Kollision stattgefunden hat,sonst true
    if(collisions.empty())
        return false;
    else
        return true;
}

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

double *RungeKuttaAdaptive::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;
}

---