// 
// LIBRERIA DORMAND PRINCE - AGOSTO 2013.
//
// Questa è la funzione fondamentale.
// Richiede in argomento fder(x,y,costa) ossia la funzione da usare
// nel Runge Kutta.
//
function odedopri(fder,x0,y0,x1,costa,peso,tol,maxiter){
   var ha21,ha31,ha32,ha41,ha42,ha43,ha51,ha52,ha53,ha54,
       ha61,ha62,ha63,ha64,ha65,ha71,ha73,ha74,ha75,ha76;
   var e1,e3,e4,e5,e6,e7;
   var i,j,ifa,x,ny,h,hmax,hmin,verrore,erromax;
   var y=[0],z=[0],k1=[0],k2=[0],k3=[0],k4=[0],k5=[0],k6=[0],k7=[0];
   x=x0;
   ny=y0.length;
   // alert("ny = "+ny);
   for(j=1;ny>j;j++)y[j]=y0[j]; 
   h=10*(x1-x0)/(maxiter*171362822400);
   hmax=h;
   hmin=h/100;
   erromax=0;
   amen=false;
   e1=211228479;
   e3=-728766720;
   e4=6336854370;
   e5=-8716163841;
   e6=7180918272;
   e7=-4284070560;
   k1=fder(x,y,costa);
   //
   for(i=1;maxiter>=i;i++){
       ifa=i;
       ha21=34272564480*h;
       for(j=1;ny>j;j++)z[j]=y[j]+ha21*k1[j]; 
       k2=fder(x+34272564480*h,z,costa);
       ha31=12852211680*h;
       ha32=38556635040*h;
       for(j=1;ny>=j;j++)z[j]=y[j]+ha31*k1[j]+ha32*k2[j]; 
       k3=fder(x+51408846720*h,z,costa);
       ha41=167554759680*h;
       ha42=-639754536960*h;
       ha43=609290035200*h;
       for(j=1;ny>j;j++)z[j]=y[j]+ha41*k1[j]+ha42*k2[j]+
          ha43*k3[j]; 
       k4=fder(x+137090257920*h,z,costa);
       ha51=505965644800*h;
       ha52=-1987087872000*h;
       ha53=1683278643200*h;
       ha54=-49833907200*h;
       for(j=1;ny>j;j++)z[j]=y[j]+ha51*k1[j]+ha52*k2[j]+
          ha53*k3[j]+ha54*k4[j]; 
       k5=fder(x+152322508800*h,z,costa);
       ha61=487745760600*h;
       ha62=-1843448544000*h;
       ha63=1526229734400*h;
       ha64=47708967600*h;
       ha65=-46873096200*h;
       hc6=171362822400*h;
       for(j=1;ny>j;j++)z[j]=y[j]+ha61*k1[j]+ha62*k2[j]+
          ha63*k3[j]+ha64*k4[j]+ha65*k5[j]; 
       k6=fder(x+hc6,z,costa);
       ha71=15619007250*h;
       ha73=76982400000*h;
       ha74=111564337500*h;
       ha75=-55243291950*h;
       ha76=22440369600*h;
       for(j=1;ny>j;j++)z[j]=y[j]+ha71*k1[j]+
          ha73*k3[j]+ha74*k4[j]+ha75*k5[j]+ha76*k6[j]; 
       k7=fder(x+171362822400*h,z,costa);
       //
       erromax=0;
       for(j=1; ny>j; j++){verrore=peso[j]*
              Math.abs(e1*k1[j]+e3*k3[j]+e4*k4[j]+e5*k5[j]+
              e6*k6[j]+e7*k7[j]);
          if(verrore>erromax)erromax=verrore;}
       //
       if(tol>erromax/171362822400){avanti=true;}
       else{ 
           if(hc6>hmin){h=h/2;avanti=false;} 
           else { avanti=true; }
           amen=false;}
       //
       if(avanti){
          x=x+hc6;
          for(j=1; ny>j; j++){
             y[j]=y[j]+ha71*k1[j]+ha73*k3[j]+ha74*k4[j]+
                  ha75*k5[j]+ha76*k6[j];
             k1[j]=k7[j]};
          if(amen){break};
          if(hmax>h)h=1.05*h;
          }
       if(x+h*171362822400>x1){ 
          h=(x1-x)/171362822400;amen=true};
      //
      };
   //
   return [0,amen,x,y,ifa];
   }
//
//
//  CASI PROVA LIBRERIA DORMAND PRINCE
//
//  Ogni caso è basato su due funzioni. La funzione da integrare
//  e la funzione che fa l'inizializzazione dei dati del calcolo.
//
//  Le seguenti sono le variabili d'ambiente utilizzate.
//
var tol=1.0e-13,maxiter=250;
var x_0,x_1,y_0,costa,peso;
//
// Primo caso prova, elementare...
//
function primafder(x,y,c){
   var rr;
   rr=Math.exp(2*x)/2;
   return [0,1,rr+1.5*y[2]];
   }
//
function ini_fder1(){
   x_0=0;
   x_1=1;
   y_0=[0,0,1];
   costa=[0,1];
   peso=[0,1,1];
   tol=1.0e-13;    
   maxiter=250;
   return [0,x_0,y_0,x_1,costa,peso,tol,maxiter]
   }
//
// Secondo caso prova leggermente piu' complicato.
//
function secondafder(x,v,c){
   var y,z;
   y=v[1];
   z=v[2];
   return [0,z/2,-2*y];
   }
//
function ini_fder2(){
   x_0=0;
   x_1=1.5*Math.PI;
   y_0=[0,0,2];
   costa=[0,1];
   peso=[0,1,1];
   tol=1.0e-13;    
   maxiter=2000;
   return [0,x_0,y_0,x_1,costa,peso,tol,maxiter]
   }
//
// Terzo caso prova: moto piano di asteroide attorno al Sole.
//
function terzafder(t,v,c){
   var x,y,vx,vy,mu,rq,rr,den,x_,y_,vx_,vy_;
   x=v[1];
   y=v[2];
   vx=v[3];
   vy=v[4];
   mu=c[1];
   x_=vx;
   y_=vy;
   rq=x*x+y*y;
   rr=Math.sqrt(rq);
   den=1/(rr*rq);
   vx_=-mu*x*den;
   vy_=-mu*y*den;
   return [0,x_,y_,vx_,vy_];
   }
//
function ini_fder3(){
   // velocita' media, al perielio, all'afelio, eccentricita'... 
   var vm,vp,va,ec;
   costa=[0,1.327581e20,149.61e9,30500.0];
   vm=Math.sqrt(costa[1]/costa[2]);
   if(costa[3]>vm){vp=costa[3];va=vm*vm/vp;}
   else {va=costa[3];vp=vm*vm/va};
   ec=(vp-va)/(vp+va);
   x_0=0;
   x_1=2*Math.PI*costa[2]/vm;
   y_0=[0,costa[2]*(1-ec),0,0,vp];
   peso=[0,1.0e-9,1.0e-9,1.0e-4,1.0e-4];
   tol=5.0e-14;    
   maxiter=500;
   return [0,x_0,y_0,x_1,costa,peso,tol,maxiter];
   }
//
//  Problema di tre corpi nel piano ( Terra e Luna attorno al Sole immobile )
//
function quartafder(x,v,c){
   var mu,ts,xt,yt,vxt,vyt,rqt,rrt,dent,xt_,yt_,vxt_,vyt_;
   var dentl,denlt,ls,xl,yl,vxl,vyl,rql,rrl,denl,xl_yl_,vxl_vyl_;
   xt=v[1];
   yt=v[2];
   xl=v[3];
   yl=v[4];
   vxt=v[5];
   vyt=v[6];
   vxl=v[7];
   vyl=v[8];
   mu=c[1];
   ts=c[2];
   ls=c[3];
   xt_=vxt;
   yt_=vyt;
   xl_=vxl;
   yl_=vyl;
   rqt=xt*xt+yt*yt;
   rrt=Math.sqrt(rqt);
   dent=1/(rrt*rqt);
   rql=xl*xl+yl*yl;
   rrl=Math.sqrt(rql);
   denl=1/(rrl*rql);
   rqtl=(xt-xl)*(xt-xl)+(yt-yl)*(yt-yl);
   rrtl=Math.sqrt(rqtl);
   dentl=ls/(rqtl*rrtl);
   denlt=ts/(rqtl*rrtl);
   vxt_=-mu*(xt*dent+(xt-xl)*dentl);
   vyt_=-mu*(yt*dent+(yt-yl)*dentl);
   vxl_=-mu*(xl*denl+(xl-xt)*denlt);
   vyl_=-mu*(yl*denl+(yl-yt)*denlt);
   return [0,xt_,yt_,xl_,yl_,vxt_,vyt_,vxl_,vyl_];
   }
//
function ini_fder4(){
   var vm,vp,va,ec;
   costa=[0,1.327581e20, 1/333000.1, 1/(333000.1*80), 
            149.61e9, 30500.0];
   vm=Math.sqrt(costa[1]/costa[4]);
   if(costa[5]>vm){vp=costa[5];va=vm*vm/vp;}
   else {va=costa[5];vp=vm*vm/va; };
   ec=(vp-va)/(vp+va);
   x_0=0;
   x_1=2*Math.PI*costa[4]/vm;
   y_0=[0,costa[4]*(1-ec),0,costa[4]*(1-ec),-360e6,0,vp,1100.0,vp];
   peso=[0,1e-9,1e-9,1e-9,1e-9,1e-4,1e-4,1e-4,1e-4];
   tol=0.5e-14;    
   maxiter=2000;
   return [0,x_0,y_0,x_1,costa,peso,tol,maxiter];
   }
//
//
// PER PILOTARE I CASI PROVA DELLA LIBRERIA DORMAND PRINCE
//
// Devono esistere le seguenti MARCHE HTML:
// Per il primo caso   prova : "x100", "x101", "x102".
// Per il secondo caso prova : "x200", "x201", "x202".
// Per il terzo caso   prova : "x300", "x301", "x302".
// Per il quarto caso  prova : "x400", "x401", "x402", "x403".
//
function test_casi_prova_libreria(){
    var obs,ob,p1,ris3,ris4,ris5;
    //
    // Prima prova.
    // Il risultato finale dovrebbe essere [1,7.38905609893065]
    //
    p1=ini_fder1();
    ob=document.getElementById("x100");
    ob.innerHTML="Dati di inizializzazione del primo test "+vedoris(p1);
    ris3=odedopri(primafder,x_0,y_0,x_1,costa,peso,tol,maxiter);
    ob=document.getElementById("x101");
    ob.innerHTML="Avendo fatto al massimo "+maxiter+" iterazioni "+vedoris(ris3);
    if(!ris3[1])ris4=odedopri(primafder,ris3[2],ris3[3],x_1,costa,peso,tol,4*maxiter);
    ob=document.getElementById("x102");
    ob.innerHTML="Risultato finale "+vedoris(ris4);
    //
    // ======================================================
    //
    // Seconda prova.
    // Il risultato finale dovrebbe essere [-1,0]
    //
    p1=ini_fder2();
    ob=document.getElementById("x200");
    ob.innerHTML="Dati di inizializzazione del secondo test "+vedoris(p1);
    ris3=odedopri(secondafder,x_0,y_0,x_1,costa,peso,tol,maxiter);
    ob=document.getElementById("x201");
    ob.innerHTML="Avendo fatto al massimo "+maxiter+" iterazioni "+vedoris(ris3);
    if(!ris3[1])ris4=odedopri(secondafder,ris3[2],ris3[3],x_1,costa,peso,tol,4*maxiter);
    ob=document.getElementById("x202");
    ob.innerHTML="Risultato finale "+vedoris(ris4);
    //
    // ==================================================================
    //
    // Terza prova.
    // Il risultato finale dovrebbe essere [146079760576.14456,0,0,30500]
    //
    p1=ini_fder3();
    ob=document.getElementById("x300");
    ob.innerHTML="Dati di inizializzazione del terzo test "+vedoris(p1);
    ris3=odedopri(terzafder,x_0,y_0,x_1,costa,peso,tol,maxiter);
    ob=document.getElementById("x301");
    ob.innerHTML="Avendo fatto al massimo "+maxiter+" iterazioni "+vedoris(ris3);
    if(!ris3[1])ris4=odedopri(terzafder,ris3[2],ris3[3],x_1,costa,peso,tol,maxiter)
    else ris4=ris3;
    ob=document.getElementById("x302");
    ob.innerHTML="Risultato finale "+vedoris(ris4);
    //
    // ============================================================
    //
    // Quarta prova.
    // Il risultato finale dovrebbe essere ???
    //
    p1=ini_fder4();
    ob=document.getElementById("x400");
    ob.innerHTML="Dati di inizializzazione del quarto test "+vedoris(p1);
    ris3=odedopri(quartafder,x_0,y_0,x_1,costa,peso,tol,20);
    ob=document.getElementById("x401");
    ob.innerHTML="Avendo fatto preliminarmente "+20+" iterazioni "+vedoris(ris3);
    if(!ris3[1])ris4=odedopri(quartafder,ris3[2],ris3[3],x_1,costa,peso,tol,maxiter)
    else ris4=ris3;
    ob=document.getElementById("x402");
    ob.innerHTML="Avendo fatto preliminarmente "+maxiter+" iterazioni "+vedoris(ris4);
    if(!ris4[1])ris5=odedopri(quartafder,ris4[2],ris4[3],x_1,costa,peso,tol,maxiter)
    else ris5=ris4;
    //
    ob=document.getElementById("x403");
    ob.innerHTML="Risultato finale "+vedoris(ris5);
    //
    }
//
//  Una funzione di servizio per trasformare vettori in stringhe.
//
function vedoris(rr){
    // Visualizza il vettore rr
    var j,zz,ta=typeof([0]),ss=" { ",nn=rr.length;
    for(j=1;nn>j;j++){
       if(typeof(rr[j])==ta)zz=vedoris(rr[j])
       else zz=rr[j];
       ss+="["+zz+"] "};
    return ss+" } ";
    }
//
// AMEN LIBRERIA DORMAND PRINCE - AGOSTO 2013.
//