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