calculator.cpp

   1 //=============================================================================
   2 //
   3 // CLASS Calculator
   4 // contains all the magic math
   5 //
   6 //=============================================================================
   7
   8
   9 //== INCLUDES =================================================================
  10
  11 #include "calculator.h"
  12
  13
  14 //== IMPLEMENTATION ===========================================================
  15
  16
  17 Calculator::Calculator(QuPlot *parent) :
  18         QDialog(parent)
  19 {
  20         c = parent;
  21         derive_init();
  22
  23 }
  24
  25 Calculator::~Calculator()
  26 {
  27 }
  28
  29
  30 //-----------------------------------------------------------------------------
  31 void Calculator::derive_init()
  32 {
  33     file_rk_bif = "../rk_bifurcate.dat";
  34     file_rk_avg = "../rk_bif_avg_velo.dat";
  35     file_rk_out = "../rk_out.dat";
  36     file_rk_pt1 = "../rk_point1.dat";
  37     file_rk_pt2 = "../rk_point2.dat";
  38     file_rk_biv = "../rk_bif_velo.dat";
  39
  40     file_eu_bif = "../eu_bifurcate.dat";
  41     file_eu_avg = "../eu_bif_avg_velo.dat";
  42     file_eu_out = "../eu_out.dat";
  43     file_eu_pt1 = "../eu_point1.dat";
  44     file_eu_pt2 = "../eu_point2.dat";
  45     file_eu_biv = "../eu_bif_velo.dat";
  46
  47     file_rk_conv = "../rk_conv.dat";
  48     file_eu_conv = "../eu_conv.dat";
  49
  50     file_rk_variance = "../rk_variance.dat";
  51     file_rk_distance = "../rk_msd.dat";
  52     v = c->valuesMap.value(1);
  53
  54 }
  55
  56 void Calculator::derivs(const Doub x, VecDoub_I &y, VecDoub_O &dydx,
  57                         const Doub params[], Normaldev &gauss)
  58 {
  59     //params[0] = a;
  60     //params[1] = A;
  61     //params[2] = F;
  62     //params[3] = l;
  63     //params[4] = w;
  64     //params[5] = eta1;
  65     //params[6] = eta2;
  66     //params[7] = l1;
  67     //params[8] = l2;
  68     //params[9] = T;
  69     //params[10] = theta;
  70     //params[11] = e0;
  71     //params[12] = e1;
  72     //params[13] = e2;
  73     //params[14] = e3;
  74     //params[15] = fx;
  75     //params[16] = fy;
  76     //static Doub G, xi, xi1, xi2, c, s, st, ax, ay, e0, e1, e2, e3, eta, l1, l2;
  77     static Doub xi, G, s, c, st, ax, ay;
  78     static Doub x1, y1, x2, y2, F1x, F1y, F2x, F2y, fx, fy, stx, sty;
  79     //for (Int i = 0;i<10;i++)
  80      //   qDebug() << "params["<<i<<"] :" << params[i];
  81
  82     switch (v) {
  83     case 1:
  84         xi = gauss.dev();
  85         G = sqrt(2.0*params[5]*params[9])*xi;
  86         dydx[0] = y[1];
  87         dydx[1] = -params[5]*y[1]-sin(y[0])+params[1]*sin(params[4]*x)+params[2]+G;
  88         break;
  89     case 2:
  90         dydx[0] = y[1];
  91         dydx[1] = -params[5]*y[1]-y[0];
  92         break;
  93     case 3:
  94         dydx[0] = y[1];
  95         dydx[1] = -y[0];
  96         break;
  97     case 4:
  98
  99         //xi1 = gauss.dev();
 100         //xi2 = gauss.dev();
 101
 102         /*
 103         eta = params[5]+params[6];
 104         e0 = 1.0/eta;
 105         e1 = 1.0/(params[5]*params[3]);
 106         e2 = 1.0/(params[6]*params[3]);
 107         e3 = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 108
 109         l1 = params[3]*e0*params[5];
 110         l2 = params[3]*e0*params[6];
 111
 112         */
 113
 114        //qDebug() << eta << "=" << params[5] << "+" << params[6];
 115        //qDebug() << e0 << "=" << "1.0/"<< eta;
 116        //qDebug() << e1 << "=" << "1.0/"<<params[5]<<"*"<<params[3];
 117        //qDebug() << e2 << "=" << "1.0/"<<params[6]<<"*"<<params[3];
 118        //qDebug() << e3 << "=" << "1.0/"<<params[5]<<"-"<<"1.0/"<<params[6];
 119
 120         //G = sqrt(2.0*eta*params[9]);
 121         G = 0;
 122
 123         c = cos(y[2]);
 124         s = sin(y[2]);
 125
 126         x1 = y[0] + params[7]*c;
 127         y1 = y[1] + params[7]*s;
 128         x2 = y[0] - params[8]*c;
 129         y2 = y[1] - params[8]*s;
 130
 131         st = sin(params[4]*x);
 132
 133         //qDebug() << st << "=" << "sin("<<params[4]<<"*"<<x<<")";
 134
 135         ax = cos(params[1])*params[0]*st;
 136         ay = sin(params[1])*params[0]*st;
 137
 138         //fx = cos(params[10])*params[2];
 139         //fy = sin(params[10])*params[2];
 140
 141         stx = ax+params[15];
 142         sty = ay+params[16];
 143
 144         //qDebug() << "ax =" << ax << "=" << params[0] << "*" << "cos("<<params[1]<<")";
 145         //qDebug() << "ay =" << ay << "=" << params[0] << "*" << "sin("<<params[1]<<")";
 146         //qDebug() << "fx =" << fx << "=" << "cos("<<params[10]<<")*"<<params[2];
 147         //qDebug() << "fy =" << fy << "=" << "sin("<<params[10]<<")*"<<params[2];
 148         //qDebug() << "stx =" << stx <<"="<< ax << "*" << st << "+"<< fx;
 149         //qDebug() << "sty =" << sty <<"="<< ay << "*" << st << "+" << fy;
 150
 151
 152
 153         //qDebug() << y[0] + params[7]*c << "=" << y[0] << "+" << params[7] << "*" << c;
 154         //qDebug() << y[1] + params[7]*c << "=" << y[1] << "+" << params[7] << "*" << s;
 155         //qDebug() << y[0] - params[8]*c << "=" << y[0] << "-" << params[8] << "*" << c;
 156         //qDebug() << y[1] - params[8]*c << "=" << y[1] << "-" << params[8] << "*" << s;
 157
 158         F1x = sin(x1)*cos(y1);//+1.6*sin(x1)+1.0*cos(y1);
 159         F1y = cos(x1)*sin(y1);//+1.6*cos(x1)+1.0*sin(y1);
 160         F2x = sin(x2)*cos(y2);//+1.6*sin(x2)+1.0*cos(y1);
 161         F2y = cos(x2)*sin(y2);//+1.6*cos(x2)+1.0*sin(y2);
 162
 163         dydx[0] = (F1x+F2x+2.0*stx)*params[11];
 164         dydx[1] = (F1y+F2y+2.0*sty)*params[11];
 165         dydx[2] = params[12]*(c*F1y-s*F1x)-params[13]*(c*F2y-s*F2x)+params[14]*(c*sty-s*stx);
 166
 167 /*
 168         c = cos(y[2]);
 169         s = sin(y[2]);
 170
 171         eta = params[5]+params[6];
 172         l1 = params[3]*params[5]/eta;
 173         l2 = params[3]*params[6]/eta;
 174
 175
 176         F1x = sin(x1)*cos(y1);//+1.6*sin(x1)+1.0*cos(y1);
 177         F1y = cos(x1)*sin(y1);//+1.6*cos(x1)+1.0*sin(y1);
 178         F2x = sin(x2)*cos(y2);//+1.6*sin(x2)+1.0*cos(y1);
 179         F2y = cos(x2)*sin(y2);//+1.6*cos(x2)+1.0*sin(y2);
 180
 181         dydx[0] = (F1x+F2x+2.0*(params[0]*cos(params[1])*sin(params[4]*x)+params[10]*params[2])/eta;
 182         dydx[1] = (F1y+F2y+2.0*(params[0]*sin(params[1])*sin(params[4]*x)+params[10]*params[2])/eta;
 183         dydx[2] = (-s*F1x+c*F1y)/(params[5]*params[3]) - (-s*F2x+c*F2y)/(params[6]*params[3])
 184                   + (1.0/params[5] - 1.0/params[6])*params[0]*( (-s*cos(params[1])+c*sin(params[1])+sin(params[4]*x))+(- );
 185 */
 186         break;
 187     }
 188 }
 189
 190
 191 //-----------------------------------------------------------------------------
 192
 193
 194 void Calculator::rk4(VecDoub_I &y, VecDoub_I &dydx, const Doub x, const Doub h,
 195                      VecDoub_O &yout, const Doub params[], Normaldev &gauss)
 196 {
 197     Int n = y.size();
 198     VecDoub dym(n),dyt(n),yt(n);
 199     Doub hh = h*0.5;
 200     Doub h6 = h/6.0;
 201     Doub xh = x+hh;
 202     for (Int i=0;i<n;i++) yt[i] = y[i]+hh*dydx[i];
 203     derivs(xh,yt,dyt,params,gauss);
 204     for (Int i=0;i<n;i++) yt[i] = y[i]+hh*dyt[i];
 205     derivs(xh,yt,dym,params,gauss);
 206     for (Int i=0;i<n;i++) {
 207         yt[i] = y[i]+h*dym[i];
 208         dym[i] += dyt[i];
 209     }
 210     derivs(x+h,yt,dyt,params,gauss);
 211     for (Int i=0;i<n;i++) {
 212         yout[i] = y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);
 213     }
 214 }
 215
 216
 217 void Calculator::euler(VecDoub_I &y, VecDoub_I &dydx, const Doub x,
 218                        const Doub h, VecDoub_O &yout)
 219 {
 220     Int n = y.size();
 221     for (Int i=0;i<n;i++) {
 222         yout[i] = y[i]+h*dydx[i];
 223     }
 224 }
 225 //-----------------------------------------------------------------------------
 226
 227
 228 /*
 229  * returns a random number between 0 and 2*pi
 230  */
 231 double Calculator::randpi(Ran &r)
 232 {
 233     Doub s = 2.0*M_PI*r.doub();
 234     return s;
 235 }
 236
 237 //-----------------------------------------------------------------------------
 238
 239 /*
 240  * Writes the of two maps into a tab-separated file
 241  * The second map *must* be specified.
 242  */
 243 void Calculator::write_to_file(const QString & name,
 244                                const QMap<double, double> & map_x,
 245                                const QMap<double, double> & map_y,
 246                                const QMap<double, double> & map_phi,
 247                                const Doub params[])
 248 {
 249     qDebug() << "writing to file" << name;
 250
 251
 252     QFile file(name);
 253     QTextStream out(&file);
 254     file.open(QIODevice::WriteOnly);
 255     out << "#";
 256     for (Int i=0;i<10;i++)
 257         out << "params["<<i<<"] = " <<  params[i] << " ";
 258     out << endl;
 259     out << "r\tx\ty\tphi" << endl;
 260     QMapIterator<double,double> x(map_x);
 261     QMapIterator<double,double> y(map_y);
 262     QMapIterator<double,double> phi(map_phi);
 263     switch(v) {
 264         case 1:
 265         while (x.hasNext()) {
 266             x.next();
 267             out << x.key() << "\t" << x.value() << endl;
 268         }
 269         break;
 270         case 4:
 271         while (x.hasNext() && y.hasNext() && phi.hasNext()) {
 272             x.next();
 273             y.next();
 274             phi.next();
 275             out << x.key() << "\t" << x.value() << "\t" << y.value() << "\t" << phi.value() << endl;
 276         }
 277         break;
 278     }
 279     file.close();
 280     qDebug () << "done.";
 281 }
 282
 283
 284 //-----------------------------------------------------------------------------
 285
 286
 287 Doub Calculator::mean(VecDoub_I &msd, const Int & samples)
 288 {
 289     Doub mu,sum,sigma;
 290     VecDoub var(samples,0.0);
 291
 292     assert(msd.size() == samples);
 293     sum = 0;
 294     for(Int i=0;i<samples;i++) {
 295         sum += msd[i];
 296     }
 297     mu = sum/samples;
 298
 299     for(Int i=0;i<samples;i++) {
 300         var[i] = SQR(msd[i])-SQR(mu);
 301     }
 302     assert(var.size() == samples) ;
 303     sum = 0;
 304     for(Int j=0;j<samples;j++) {
 305         sum += var[j];
 306     }
 307
 308     sigma = sum/samples;
 309
 310     //qDebug() << sigma;
 311     return sigma;
 312 }
 313
 314 /*
 315  * This function is kinda redundant right now, but will be needed later
 316  */
 317 void Calculator::gen_gnuplot(const QString & name)
 318 {
 319
 320     qDebug() << "generating gnuplot file";
 321
 322     QFile file(name);
 323     QTextStream out(&file);
 324     file.open(QIODevice::WriteOnly);
 325
 326     out << "L = 2*pi;";
 327     out << "unset key;";
 328     out << "set lmargin 10;";
 329     out << "set rmargin 3;";
 330     out << "set multiplot;";
 331     out << "set tmargin = 3;";
 332     out << "set bmargin = 0;";
 333     out << "set border 2+4+8;";
 334     out << "set size 1,0.66667;";
 335     out << "set origin 0.0,0.33333;";
 336     out << "set xtics;";
 337     out << "set format x \"\";";
 338     out << "unset xlabel;";
 339     out << "set ylabel \"x(kT) mod L \";";
 340     out << "set ytics (\"0\" 0, \"L/2\", \"L\" L);";
 341     out << "set yrange [0:L];";
 342     out << "plot \""<<file_rk_bif<<"\" w dots 1;";
 343     out << "unset title;";
 344     out << "set tmargin 0;";
 345     out << "set bmargin 3;";
 346     out << "unset format;";
 347     out << "set size 1,0.33333;";
 348     out << "set border 1+2+4+8;";
 349     out << "set origin 0.0,0.0;";
 350     out << "set xtics 0,0.1;";
 351     out << "set xlabel \"a\";";
 352     out << "set ylabel \"v\";";
 353     out << "set yrange [-0.75:0.75];";
 354     out << "set ytics -0.5,0.5;";
 355     out << "plot \""<<file_rk_avg<<"\" w dots 3;";
 356     out << "unset multiplot;";
 357
 358 }
 359
 360 void Calculator::filter(QMap<double, double> & map, const Doub & eps)
 361 {
 362     qDebug() << "filter crap smaller than " << eps;
 363     Int size, k;
 364     Doub diff, key;
 365     QList<double> values, list;
 366     list = map.uniqueKeys();
 367     //qDebug() << list;
 368     size = list.size();
 369     //qDebug() << "size of map:" << size;
 370     for(Int i = 0; i < size; i++) {
 371         key = list.at(i);
 372         //qDebug() << "["<<key<<"]" << "writing values into list";
 373         values = map.values(key);
 374         //qDebug() << "["<<key<<"]" << values.size() << "values for" << key;
 375         k = 0;
 376         for(Int k=0;k<values.size()-1;k++) {
 377             diff = abs(values.at(k+1) - values.at(k));
 378             //qDebug() << values.at(k+1) <<"-"<< values.at(k) << "=" << diff;
 379             if (diff <= eps) {
 380                 //qDebug() << diff << "<=" << eps;
 381                 //qDebug() << "["<<k<<"]" << "removing" << values.at(k) << "from list";
 382                 values.replace(k+1,values.at(k));
 383             }
 384         } // end for
 385         for (int j = 0; j < values.size(); j++)
 386             //qDebug() << "["<<j<<"]" << values.at(j);
 387         //qDebug() << "removing" << key << "from map";
 388         map.remove(key);
 389         for (int j = 0; j < values.size(); j++) {
 390             //qDebug() << "inserting" << values.at(j) << "into map at" << key;
 391             map.insertMulti(key,values.at(j));
 392         }
 393         //qDebug() << map.values(key);
 394         //qDebug() << endl << "end while" << endl;
 395     }
 396 }
 397
 398 void Calculator::IC_setup(VecDoub &y, Doub x0, Doub xdot0, Doub y0, Doub ydot0,
 399                           Doub phi0, Doub params[])
 400 {
 401     Int n = y.size();
 402     const Int sys = c->valuesMap.value(1);
 403
 404     /*
 405     qDebug() << "using " << "\tx0 = " << x0 << "\ty0 = " << y0
 406             << "\txdot0 = " << xdot0 << "\tydot0 = " << ydot0
 407             << "\tphi0 = " << fmod(phi0,360) << endl;
 408     */
 409
 410     switch(sys) {
 411         case 1:
 412         params[0] = c->valuesMap.value(2); // a;
 413         params[1] = c->valuesMap.value(3); // A;
 414         params[2] = c->valuesMap.value(4); // F;
 415         params[3] = c->valuesMap.value(5); // l;
 416         params[4] = c->valuesMap.value(6); // w;
 417         params[5] = c->valuesMap.value(25); // eta1;
 418         params[6] = 0;
 419         params[9] = c->valuesMap.value(32); // T;
 420         params[10] = 0;
 421         y[0] = x0, y[1] = xdot0;
 422         break;
 423         case 4:
 424         params[0] = c->valuesMap.value(2); // a;
 425         params[1] = fmod(c->valuesMap.value(3),360)*PI/180; // alpha;
 426         params[2] = c->valuesMap.value(4); // F;
 427         params[3] = c->valuesMap.value(5); // l;
 428         params[4] = c->valuesMap.value(6); // w;
 429         params[5] = c->valuesMap.value(25); // eta1;
 430         params[6] = c->valuesMap.value(26); // eta2;
 431         params[9] = c->valuesMap.value(32); // T;
 432         params[10] = fmod(c->valuesMap.value(40),360)*PI/180; // theta;
 433         Doub e0 = 1.0/(params[5]+params[6]);
 434         Doub e1 = 1.0/(params[5]*params[3]);
 435         Doub e2 = 1.0/(params[6]*params[3]);
 436         Doub e3 = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 437         Doub fx = cos(params[10])*params[2];
 438         Doub fy = sin(params[10])*params[2];
 439         params[11] = e0; // 1/eta
 440         params[12] = e1;
 441         params[13] = e2;
 442         params[14] = e3;
 443         params[15] = fx;
 444         params[16] = fy;
 445         params[7] = params[3]*params[11]*params[5]; // l1;
 446         params[8] = params[3]*params[11]*params[6]; // l2;
 447         y[0] = x0, y[1] = y0, y[2] = phi0;
 448         break;
 449     }
 450     /*
 451     for (Int i = 0; i<n;i++)
 452         qDebug() << "y("<<i<<") =" << y[i];
 453     qDebug() << "eta1 = " << eta1 << "\teta2 = " << eta2 << "\tA = " << A << "\tF = " << F;
 454     */
 455 }
 456
 457 void Calculator::bifparam(const int &v, const QString &param, const Doub &r,
 458                           Doub params[])
 459 {
 460     switch(v) {
 461         case 1:
 462         if (param == "E")
 463             params[5] = r;
 464         else if (param == "F")
 465             params[2] = r;
 466         else if (param == "A")
 467             params[1] = r;
 468         else
 469             throw ("Fatal error!");
 470         break;
 471         case 4:
 472         if (param == "A")
 473             params[1] = fmod(r,360)*PI/180;
 474         else if (param == "F") {
 475             params[2] = r;
 476             params[15] = cos(params[10])*params[2];
 477             params[16] = sin(params[10])*params[2];
 478         }
 479         else if (param == "a")
 480             params[0] = r;
 481         else if (param == "w")
 482             params[4] = r;
 483         else if (param == "l") {
 484             params[3] = r;
 485             params[11] = 1.0/(params[5]+params[6]);
 486             params[12] = 1.0/(params[5]*params[3]);
 487             params[13] = 1.0/(params[6]*params[3]);
 488             params[14] = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 489             params[7] = params[3]*params[11]*params[5]; // l1;
 490             params[8] = params[3]*params[11]*params[6]; // l2;
 491         }
 492         else if (param == "e") {
 493             params[5] = r;
 494             params[11] = 1.0/(params[5]+params[6]);
 495             params[12] = 1.0/(params[5]*params[3]);
 496             params[13] = 1.0/(params[6]*params[3]);
 497             params[14] = 1.0/params[3]*((1.0/params[5])-(1.0/params[6]));
 498             params[7] = params[3]*params[11]*params[5]; // l1;
 499             params[8] = params[3]*params[11]*params[6]; // l2;
 500         }
 501         else if (param == "t")
 502             params[10] = fmod(r,360)*PI/180;
 503         else
 504             throw("Fatal error!");
 505         break;
 506         qDebug() << "Bifurcation parameter is " << param;
 507     }
 508 }
 509
 510 void Calculator::set_color(const Doub vel, QString & color)
 511 {
 512     Doub ctol = 0.05;
 513     if (vel > 1.0-ctol && vel < 1.0+ctol)
 514         color = "red";
 515     else if (vel > 0.5-ctol && vel < 0.5+ctol)
 516         color = "orange";
 517     else if (vel > 0.0-ctol && vel < 0.0+ctol)
 518         color = "grey";
 519     else if (vel > -0.5-ctol && vel < -0.5+ctol)
 520         color = "green";
 521     else if (vel > -1.0-ctol && vel < -1.0+ctol)
 522         color = "blue";
 523     else color = "black";
 524 }
 525
 526 void Calculator::calc2(const Int n)
 527 {
 528     // setting up a few start values
 529     qDebug() << "calc2 has been called: " << "calc2("<<n<<")";
 530     Int per;
 531     Normaldev gauss(0,1,time(NULL));
 532     Ran ran(time(NULL));
 533     Doub x1,x2,y1,y2,xtest;
 534     Doub L, tau;
 535     Doub x0, xdot0, y0, ydot0, phi0;
 536     x0 = c->valuesMap.value(8);
 537     xdot0 = c->valuesMap.value(10);
 538     ydot0 = c->valuesMap.value(24);
 539     y0 = c->valuesMap.value(23);
 540     phi0 = c->valuesMap.value(27);
 541     const Int samples = c->valuesMap.value(35);
 542     const Int v = c->valuesMap.value(1);
 543     const Int M = c->valuesMap.value(11);
 544     const Doub fromt = c->valuesMap.value(22);
 545     const Int s = c->valuesMap.value(14);
 546     const Doub h = c->valuesMap.value(13);
 547     const Doub t0 = c->valuesMap.value(7);
 548     const Doub fromr1 = c->valuesMap.value(17);
 549     const Doub fromr2 = c->valuesMap.value(37);
 550     const Doub tor1 = c->valuesMap.value(18);
 551     const Doub tor2 = c->valuesMap.value(38);
 552     const Doub stepr1 = c->valuesMap.value(19);
 553     const Doub stepr2 = c->valuesMap.value(39);
 554     const Doub eps = c->valuesMap.value(36);
 555     const Doub from = c->valuesMap.value(15);
 556     const Doub to = c->valuesMap.value(16);
 557     Doub params[256];
 558     //const Doub L = 1;
 559     const bool ruk = c->ruk;
 560     const bool eul = c->eul;
 561     QString param1 = c->str1;
 562     QString param2 = c->str2;
 563     qDebug() << "System" << v;
 564
 565     Doub x;
 566
 567     Doub xpi[samples][s];
 568     Doub ypi[samples][s];
 569     Doub phipi[samples][s];
 570     VecDoub y(n,0.0);
 571     VecDoub dydx(n,0.0);
 572     VecDoub yout(n,0.0);
 573     VecDoub avg(n,0.0);
 574
 575     if (ruk) {
 576         switch(M) {
 577             case 0:
 578             qDebug() << "standard mode";
 579             qDebug() << "stepsize:" << h;
 580             qDebug() << "number of steps:" << s;
 581             x = t0;
 582             IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 583             per = 0;
 584             L = 2*PI;
 585             tau = 2*PI/params[4];
 586             while(per < s) {
 587                 xtest = x+h;
 588                 if (xtest >= per*tau) {
 589                     per++;
 590                 }
 591                 derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4()
 592                 rk4(y,dydx,x,h,yout,params,gauss);
 593                 if (per >= s*fromt) {
 594                     switch(v) {
 595                         case 1:
 596                         rkValues_x.insert(x,y[0]);
 597                         break;
 598                         case 4:
 599                         //xpi = ABS(fmod(yout[0],(2*PI)));
 600                         //ypi = ABS(fmod(yout[1],(2*PI)));
 601                         rkValues_x.insert(x,y[0]);
 602                         rkValues_y.insert(x,y[1]);
 603                         //rkValues_x.insert(x,xpi);
 604                         //rkValues_y.insert(x,ypi);
 605                         rkValues_phi.insert(x,y[2]);
 606                         x1 = y[0] + params[7]*cos(y[2]);
 607                         y1 = y[1] + params[7]*sin(y[2]);
 608                         x2 = y[0] - params[8]*cos(y[2]);
 609                         y2 = y[1] - params[8]*sin(y[2]);
 610                         rkValues_p1x.insert(x,x1);
 611                         rkValues_p1y.insert(x,y1);
 612                         rkValues_p2x.insert(x,x2);
 613                         rkValues_p2y.insert(x,y2);
 614                         break;
 615                     }
 616                 }
 617                 for (Int j=0;j<n;j++)
 618                     y[j] = yout[j];
 619                 x += h;
 620             } // end while
 621             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 622             qDebug() << "sizeof(rkValues_y) =" << rkValues_y.size();
 623             qDebug() << "sizeof(rkValues_phi) =" << rkValues_phi.size();
 624             write_to_file(file_rk_out, rkValues_x, rkValues_y, rkValues_phi, params);
 625             write_to_file(file_rk_pt1, rkValues_p1x, rkValues_p1y, rkValues_phi, params);
 626             write_to_file(file_rk_pt2, rkValues_p2x, rkValues_p2y, rkValues_phi, params);
 627             break;
 628             case 1:
 629             qDebug() << "converge mode";
 630             Int iter;
 631             Doub hh, ss;
 632             IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 633             for (iter=1;(iter*from)<to;iter++) {
 634                 hh = (iter*from);
 635                 ss = (s/hh);
 636                 x = t0;
 637                 qDebug() << "stepsize is" << hh << "with" << ss << "total steps.";
 638                 IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 639                 for (Int i=0;i<ss&&x<(ss-0.5);i++) {
 640                     derivs(x,y,dydx,params,gauss);
 641                     rk4(y,dydx,x,hh,yout,params,gauss);
 642                     for (Int j=0;j<n;j++)
 643                         y[j] = yout[j];
 644                     x += hh;
 645                 }
 646                 switch(v) {
 647                     case 1:
 648                     rkValues_x.insert(hh,yout[0]);
 649                     break;
 650                     case 4:
 651                     rkValues_x.insert(hh,yout[0]);
 652                     rkValues_y.insert(hh,yout[1]);
 653                     rkValues_phi.insert(hh,yout[2]);
 654                     break;
 655                 }
 656                 x = t0;
 657                 IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 658                 for (Int i=0;i<ss&&x<(ss-0.5);i++) {
 659                     derivs(x,y,dydx,params,gauss);
 660                     euler(y,dydx,x,hh,yout);
 661                     for (Int j=0;j<n;j++)
 662                         y[j] = yout[j];
 663                     x += hh;
 664                 }
 665                 switch(v) {
 666                     case 1:
 667                     euValues_x.insert(hh,yout[0]);
 668                     break;
 669                     case 4:
 670                     euValues_x.insert(hh,yout[0]);
 671                     euValues_y.insert(hh,yout[1]);
 672                     euValues_phi.insert(hh,yout[2]);
 673                     break;
 674                 }
 675             }
 676             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 677             write_to_file(file_rk_conv, rkValues_x, rkValues_y, rkValues_phi, params);
 678             write_to_file(file_eu_conv, euValues_x, euValues_y, euValues_phi, params);
 679             break;
 680             case 2:
 681             qDebug() << "bifurcation mode";
 682             qDebug() << "bifurcation parameter: " << param1;
 683             qDebug() << "Number of samples : " << samples;
 684             Doub r;
 685             for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
 686                 r = fromr1+j*stepr1;
 687                 qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
 688                 for (Int z=0;z<samples;z++) {
 689                     x = t0;
 690                     x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 691                     //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
 692                     IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 693                     bifparam(v,param1,r,params);
 694                     //qDebug() << x << y[0] << y[1] << y[2];
 695                     //for(Int i = 0;i<=10;i++)
 696                      //   qDebug() << params[i];
 697                     // inner for: iterate over runge-kutta
 698                     Int k = 0;
 699                     Doub hnew;
 700                     L = 2*PI;
 701                     tau = 2*PI/params[4];
 702                     //qDebug() << "tau: " << tau << "\tL: " << L;
 703                     while (k<s) {
 704                         xtest = x+h; // test next x increment
 705                         //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 706                         if (xtest >= k*tau) {
 707                             //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 708                             //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 709                             hnew = k*tau-x;
 710                             //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 711                             derivs(x,y,dydx,params,gauss);
 712                             rk4(y,dydx,x,hnew,yout,params,gauss);
 713                             switch(v) {
 714                                 case 1:
 715                                 xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 716                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 717                                 if (k >= s*fromt) {
 718                                     rkValues_x.insert(r,xpi[z][k]);
 719                                     rkValues_xdot.insert(r,dydx[0]);
 720                                     }
 721                                 break;
 722                                 case 4:
 723                                 xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 724                                 //qDebug() << "xpi["<<z<<"]["<<k<<"]:" << xpi[z][k];
 725                                 ypi[z][k] = ABS(fmod(yout[1],(2*PI)));
 726                                 phipi[z][k] = ABS(fmod(yout[2],(2*PI)));
 727                                 rkDist_x.insert(r,SQR(yout[0]-x0));
 728                                 rkDist_y.insert(r,SQR(yout[1]-y0));
 729                                 rkDist_phi.insert(r,SQR(yout[2]-phi0));
 730                                 if (k >= s*fromt) {
 731                                     rkValues_x.insert(r,xpi[z][k]);
 732                                     rkValues_y.insert(r,ypi[z][k]);
 733                                     rkValues_phi.insert(r,yout[2]);
 734                                     rkValues_xdot.insert(r,dydx[0]);
 735                                     rkValues_ydot.insert(r,dydx[1]);
 736                                     rkValues_phidot.insert(r,dydx[2]);
 737                                 }
 738                                 break;
 739                             } // end switch
 740                             k++;
 741                         }
 742                         derivs(x,y,dydx,params,gauss);
 743                         rk4(y,dydx,x,h,yout,params,gauss); // xold
 744                         for (Int j=0;j<n;j++)
 745                             y[j] = yout[j]; // y[x+h] !!!!!!!
 746                         x += h; // do it for real
 747                     } // end while
 748                     //qDebug() << "at" << x << "in time for k = " << k;
 749                     switch(v) {
 750                         Int m_x, m_y, m_phi;
 751                         case 1:
 752                         m_x = yout[0]/L;
 753                         avg[0] = (Doub) m_x/k;
 754                         rkAvgvelo_x.insert(r,avg[0]);
 755                         break;
 756                         case 4:
 757                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
 758                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
 759                         rkAvgvelo_x.insert(r,avg[0]);
 760                         rkAvgvelo_y.insert(r,avg[1]);
 761                         rkAvgvelo_phi.insert(r,avg[2]);
 762                         break;
 763                     } // end switch
 764                 } // samples loop
 765                 VecDoub msd_x(samples,0.0);
 766                 VecDoub msd_y(samples,0.0);
 767                 VecDoub msd_phi(samples,0.0);
 768                 VecDoub sigma_x(s,0.0);
 769                 VecDoub sigma_y(s,0.0);
 770                 VecDoub sigma_phi(s,0.0);
 771                 //qDebug() << "samples:" << samples << "steps:" << s;
 772                 for(Int i=0;i<s;i++) {
 773                     //qDebug() << "["<<i<<"]:-----------------------------------------";
 774                     for(Int j=0;j<samples;j++) {
 775                         //qDebug() << "xpi["<<j<<"]["<<i<<"]:" << xpi[j][i];
 776                         msd_x[j] = xpi[j][i];
 777                         //qDebug() << "["<<j<<"]:" << sper_x[j];
 778                         msd_y[j] = ypi[j][i];
 779                         msd_phi[j] = phipi[j][i];
 780                     }
 781                     sigma_x[i] = mean(msd_x,samples);
 782                     sigma_y[i] = mean(msd_y,samples);
 783                     sigma_phi[i] = mean(msd_phi,samples);
 784                     //qDebug() << "sigma_x["<<i<<"]:" << sigma_x[i];
 785                 }
 786                 for(Int k=0;k<s;k++) {
 787                     rkVariance_x.insert(r,sigma_x[k]);
 788                     rkVariance_y.insert(r,sigma_y[k]);
 789                     rkVariance_phi.insert(r,sigma_phi[k]);
 790                 }
 791                //qDebug() << rkVariance_x;
 792             } // end outer for
 793             /*
 794              * we need to do something about this
 795              */
 796             /*
 797             filter(rkValues_x,eps);
 798             filter(rkValues_y,eps);
 799             filter(rkValues_xdot,eps);
 800             filter(rkValues_ydot,eps);
 801             filter(rkAvgvelo_x, eps);
 802             filter(rkAvgvelo_y, eps);
 803             filter(rkAvgvelo_phi,eps);*/
 804             qDebug() << "sizeof(rkValues_x) =" << rkValues_x.size();
 805             qDebug() << "sizeof(rkValues_y) =" << rkValues_y.size();
 806             qDebug() << "sizeof(rkValues_xdot) =" << rkValues_xdot.size();
 807             qDebug() << "sizeof(rkValues_ydot) =" << rkValues_ydot.size();
 808             write_to_file(file_rk_bif, rkValues_x, rkValues_y, rkValues_phi, params);
 809             write_to_file(file_rk_biv, rkValues_xdot, rkValues_ydot, rkValues_phidot, params);
 810             write_to_file(file_rk_avg, rkAvgvelo_x, rkAvgvelo_y, rkAvgvelo_phi, params);
 811             write_to_file(file_rk_variance, rkVariance_x, rkVariance_y, rkVariance_phi,params);
 812             write_to_file(file_rk_distance, rkDist_x, rkDist_y, rkDist_phi, params);
 813             break;
 814             case 3:
 815             qDebug() << "dual bifurcation mode";
 816             qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
 817             qDebug() << "Number of samples : " << samples;
 818             Doub stp1, stp2;
 819             QString vx_color, vy_color;
 820             QString name = "dual_bif.dat";
 821             QFile file(name);
 822             QTextStream out(&file);
 823             file.open(QIODevice::WriteOnly);
 824             out << "x\ty\tvx\tvx_color\tvy\tvy_color\tvphi" << endl;
 825             for(Int it1=0;(it1*stepr1)<=(tor1-fromr1);it1++) {
 826                 for(Int it2=0;(it2*stepr2)<=(tor2-fromr2);it2++) {
 827                     stp1 = fromr1+it1*stepr1;
 828                     stp2 = fromr2+it2*stepr2;
 829                     qDebug() << stp1 << "("<<tor1<<")" << stp2 << "("<<tor2<<")";
 830                     for (Int z=0;z<samples;z++) {
 831                         x = t0;
 832                         x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 833                         IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 834                         bifparam(v,param1,stp1,params);
 835                         bifparam(v,param2,stp2,params);
 836                         /*
 837                         out << "# ";
 838                         for (Int i=0;i<10;i++)
 839                             out << "params["<<i<<"] = " <<  params[i] << " ";
 840                         out << endl;
 841                         */
 842                         // inner for: iterate over runge-kutta
 843                         Int k = 0;
 844                         Doub hnew;
 845                         L = 2*PI;
 846                         tau = 2*PI/params[4];
 847                         //qDebug() << "tau: " << tau << "\tL: " << L;
 848                         while (k<s) {
 849                             xtest = x+h; // test next x increment
 850                             //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 851                             if (xtest >= k*tau) {
 852                                 //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 853                                 //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 854                                 hnew = k*tau-x;
 855                                 //qDebug() << "making rk4() with stepsize"<<hnew<<"at"<<r;
 856                                 derivs(x,y,dydx,params,gauss);
 857                                 rk4(y,dydx,x,hnew,yout,params,gauss);
 858                                 k++;
 859                             }
 860                             derivs(x,y,dydx,params,gauss);
 861                             rk4(y,dydx,x,h,yout,params,gauss); // xold
 862                             for (Int j=0;j<n;j++)
 863                                 y[j] = yout[j]; // y[x+h] !!!!!!!
 864                             x += h; // do it for real
 865                         } // end while
 866                         Int m_x, m_y, m_phi;
 867                         m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
 868                         avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k;
 869                         avg[2] = (Doub) m_phi/k;
 870                         set_color(avg[0], vx_color);
 871                         set_color(avg[1], vy_color);
 872                         out << stp1 << "\t" << stp2 << "\t"
 873                                 << avg[0] <<  "\t" << vx_color << "\t" << avg[1]
 874                                 << "\t" << vy_color << "\t" << avg[2] << endl;
 875                     } // end samples
 876                 } // end 2. outer for
 877             } // end 1. outer for
 878             file.close();
 879             break;
 880         } // end switch
 881     } // end if
 882
 883     if (eul) {
 884         switch(M) {
 885           case 0:
 886           qDebug() << "(euler) standard mode";
 887           qDebug() << "stepsize:" << h;
 888           qDebug() << "number of steps:" << s;
 889           x = t0;
 890           IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 891           per = 0;
 892           L = 2*PI;
 893           tau = 2*PI/params[4];
 894           while(per < s) {
 895               xtest = x+h;
 896               if (xtest >= per*tau) {
 897                   per++;
 898               }
 899               derivs(x,y,dydx,params,gauss); // sets first dydx[] for rk4()
 900               euler(y,dydx,x,h,yout);
 901               if (per >= s*fromt) {
 902                   switch(v) {
 903                       case 1:
 904                       euValues_x.insert(x,y[0]);
 905                       break;
 906                       case 4:
 907                       //xpi = ABS(fmod(yout[0],(2*PI)));
 908                       //ypi = ABS(fmod(yout[1],(2*PI)));
 909                       euValues_x.insert(x,y[0]);
 910                       euValues_y.insert(x,y[1]);
 911                       //rkValues_x.insert(x,xpi);
 912                       //rkValues_y.insert(x,ypi);
 913                       euValues_phi.insert(x,y[2]);
 914                       x1 = y[0] + params[7]*cos(y[2]);
 915                       y1 = y[1] + params[7]*sin(y[2]);
 916                       x2 = y[0] - params[8]*cos(y[2]);
 917                       y2 = y[1] - params[8]*sin(y[2]);
 918                       euValues_p1x.insert(x,x1);
 919                       euValues_p1y.insert(x,y1);
 920                       euValues_p2x.insert(x,x2);
 921                       euValues_p2y.insert(x,y2);
 922                       break;
 923                   }
 924               }
 925               for (Int j=0;j<n;j++)
 926                   y[j] = yout[j];
 927               x += h;
 928           } // end while
 929           qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
 930           qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
 931           qDebug() << "sizeof(euValues_phi) =" << euValues_phi.size();
 932           write_to_file(file_eu_out, euValues_x, euValues_y, euValues_phi, params);
 933           write_to_file(file_eu_pt1, euValues_p1x, euValues_p1y, euValues_phi, params);
 934           write_to_file(file_eu_pt2, euValues_p2x, euValues_p2y, euValues_phi, params);
 935           break;
 936           case 2:
 937           qDebug() << "bifurcation mode";
 938           qDebug() << "bifurcation parameter: " << param1;
 939           qDebug() << "Number of samples : " << samples;
 940           Doub r;
 941           for (Int j=0;(j*stepr1)<=(tor1-fromr1);j++) {
 942               r = fromr1+j*stepr1;
 943               qDebug() << "at step" << j << "of" << (tor1-fromr1)/stepr1;
 944               for (Int z=0;z<samples;z++) {
 945                   x = t0;
 946                   x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
 947                   //qDebug() << "sample" << z << "using" << x0 << xdot0 << y0 << ydot0 << phi0;
 948                   IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
 949                   bifparam(v,param1,r,params);
 950                   // inner for: iterate over runge-kutta
 951                   Int k = 0;
 952                   Doub hnew;
 953                   L = 2*PI;
 954                   tau = 2*PI/params[4];
 955                   //qDebug() << "tau: " << tau << "\tL: " << L;
 956                   while (k<s) {
 957                       xtest = x+h; // test next x increment
 958                       //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
 959                       if (xtest >= k*tau) {
 960                           //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
 961                           //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
 962                           hnew = k*tau-x;
 963                           //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
 964                           derivs(x,y,dydx,params,gauss);
 965                           euler(y,dydx,x,hnew,yout);
 966                           if (k >= s*fromt) {
 967                               switch(v) {
 968                                   case 1:
 969                                   xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 970                                   euValues_x.insert(r,xpi[z][k]);
 971                                   euValues_xdot.insert(r,dydx[0]);
 972                                   break;
 973                                   case 4:
 974                                   xpi[z][k] = ABS(fmod(yout[0],(2*PI)));
 975                                   ypi[z][k] = ABS(fmod(yout[1],(2*PI)));
 976                                   euValues_x.insert(r,xpi[z][k]);
 977                                   euValues_y.insert(r,ypi[z][k]);
 978                                   euValues_phi.insert(r,yout[2]);
 979                                   euValues_xdot.insert(r,dydx[0]);
 980                                   euValues_ydot.insert(r,dydx[1]);
 981                                   euValues_phidot.insert(r,dydx[2]);
 982                                   break;
 983                               }
 984                           }
 985                           k++;
 986                       }
 987                       derivs(x,y,dydx,params,gauss);
 988                       euler(y,dydx,x,h,yout); // xold
 989                       for (Int j=0;j<n;j++)
 990                           y[j] = yout[j]; // y[x+h] !!!!!!!
 991                       x += h; // do it for real
 992                   } // end while
 993                   //qDebug() << "at" << x << "in time for k = " << k;
 994                   switch(v) {
 995                       Int m_x, m_y, m_phi;
 996                       case 1:
 997                       m_x = yout[0]/L;
 998                       avg[0] = (Doub) m_x/k;
 999                       euAvgvelo_x.insert(r,avg[0]);
1000                       break;
1001                       case 4:
1002                       m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
1003                       avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
1004                       euAvgvelo_x.insert(r,avg[0]);
1005                       euAvgvelo_y.insert(r,avg[1]);
1006                       euAvgvelo_phi.insert(r,avg[2]);
1007                       break;
1008                   } // end switch
1009               } // samples loop
1010           } // end outer for
1011           /*
1012            * we need to do something about this
1013            */
1014           filter(euValues_x,eps);
1015           filter(euValues_y,eps);
1016           filter(euValues_xdot,eps);
1017           filter(euValues_ydot,eps);
1018           filter(euAvgvelo_x, eps);
1019           filter(euAvgvelo_y, eps);
1020           filter(euAvgvelo_phi,eps);
1021           qDebug() << "sizeof(euValues_x) =" << euValues_x.size();
1022           qDebug() << "sizeof(euValues_y) =" << euValues_y.size();
1023           qDebug() << "sizeof(euValues_xdot) =" << euValues_xdot.size();
1024           qDebug() << "sizeof(euValues_ydot) =" << euValues_ydot.size();
1025           write_to_file(file_eu_bif, euValues_x, euValues_y, euValues_phi, params);
1026           write_to_file(file_eu_biv, euValues_xdot, euValues_ydot, euValues_phidot, params);
1027           write_to_file(file_eu_avg, euAvgvelo_x, euAvgvelo_y, euAvgvelo_phi, params);
1028           break;
1029           case 3:
1030           qDebug() << "dual bifurcation mode";
1031           qDebug() << "bifurcation parameters: " << param1<<"\t"<<param2;
1032           qDebug() << "Number of samples : " << samples;
1033           Doub stp1,stp2;
1034           QString name = "dual_bif.dat";
1035           QFile file(name);
1036           QTextStream out(&file);
1037           file.open(QIODevice::WriteOnly);
1038           for (Int it1=0;(it1*stepr1)<=(tor1-fromr1);it1++) {
1039               for (Int it2=0;(it2*stepr2)<=(tor2-fromr1);it2++) {
1040                   stp1 = fromr1+it1*stepr1;
1041                   stp2 = fromr2+it2*stepr2;
1042                   qDebug() << stp1 << "("<<tor1<<")" << stp2 << "("<<tor2<<")";
1043                   for (Int z=0;z<samples;z++) {
1044                       x = t0;
1045                       x0 = randpi(ran), xdot0 = randpi(ran), y0 = randpi(ran), ydot0 = randpi(ran), phi0 = randpi(ran);
1046                       IC_setup(y,x0,xdot0,y0,ydot0,phi0,params);
1047                       bifparam(v,param1,stp1,params);
1048                       bifparam(v,param2,stp2,params);
1049                       // inner for: iterate over runge-kutta
1050                       Int k = 0;
1051                       Doub hnew;
1052                       L = 2*PI;
1053                       tau = 2*PI/params[4];
1054                       //qDebug() << "tau: " << tau << "\tL: " << L;
1055                       while (k<s) {
1056                           xtest = x+h; // test next x increment
1057                           //if (xtest < k*tau) { qDebug() << "["<<k<<"]" << xtest << "<" << k*tau << "..."; }
1058                           if (xtest >= k*tau) {
1059                               //qDebug() << "["<<k<<"]" << xtest << ">" << k*tau << "!";
1060                               //qDebug() << "new h is" << k << "*" << tau << "-" << x << "=" << hnew;
1061                               hnew = k*tau-x;
1062                               //qDebug() << "making euler() with stepsize"<<hnew<<"at"<<r;
1063                               derivs(x,y,dydx,params,gauss);
1064                               euler(y,dydx,x,hnew,yout);
1065                               k++;
1066                           }
1067                           derivs(x,y,dydx,params,gauss);
1068                           euler(y,dydx,x,h,yout); // xold
1069                           for (Int j=0;j<n;j++)
1070                               y[j] = yout[j]; // y[x+h] !!!!!!!
1071                           x += h; // do it for real
1072                       } // end while
1073                       Int m_x, m_y, m_phi;
1074                       m_x = yout[0]/L; m_y = yout[1]/L; m_phi = yout[2]/L;
1075                       avg[0] = (Doub) m_x/k; avg[1] = (Doub) m_y/k; avg[2] = (Doub) m_phi/k;
1076                       out << scientific << stp1 << "\t" << stp2 << "\t" << avg[0] << "\t" << avg[1] << "\t" << avg[2] << endl;
1077                   } // end samples
1078               } // end 2. outer for
1079           } // end 1. outer for
1080           file.close();
1081           break;
1082         } // end switch
1083     } // end if
1084
1085     if (!ruk && !eul)
1086         throw("No method specified!");
1087
1088     rkValues_x.clear();
1089     rkValues_xdot.clear();
1090     rkValues_y.clear();
1091     rkValues_ydot.clear();
1092     rkValues_phi.clear();
1093     rkValues_phidot.clear();
1094     rkAvgvelo_x.clear();
1095     rkAvgvelo_y.clear();
1096     rkAvgvelo_phi.clear();
1097
1098     euValues_x.clear();
1099     euValues_xdot.clear();
1100     euValues_y.clear();
1101     euValues_ydot.clear();
1102     euValues_phi.clear();
1103     euValues_phidot.clear();
1104     euAvgvelo_x.clear();
1105     euAvgvelo_y.clear();
1106     euAvgvelo_phi.clear();
1107
1108 } // end calculator
1109
1110 //=============================================================================