//
//   Libreria ECMAscript del Runge Kutta di J.R.Cash e A.H.Karp.
//   Versione maggio 2006.
//
function cki(h,limite){
     var j,delta,erromax,n=y0.length;
     //
     // h rappresenta il passo diviso 979292160 ovvero 4096*27*5*7*11*23.
     // limite rappresenta il massimo errore stimato consentito.
     // I vettori y0,yn,k1,k2,k3,k4,k5,k6 debbono essere definiti 
     // esternamente e debbono essere tutti di uguale lunghezza.
     // La funzione fader(y) riceve in argomento il vettore delle
     // variabili e produce come risultato il vettore delle derivate prime.
     // Il tempo non è trattato separatamente ma va eventualmente incluso
     // come variabile la cui derivata prima è la costante 1.
     //
     // Se il risultato è false bisogna ridurre h ossia il passo. 
     //
     var a1=195858432*h;
     var b1=73446912*h,b2=220340736*h;
     var c1=293787648*h,c2=-881362944*h,c3=1175150592*h;
     var d1=-199485440*h,d2=2448230400*h,d3=-2538905600*h,
         d4=1269452800*h;
     var e1=28885010*h,e2=334719000*h,e3=40733000*h,
         e4=392055125*h,e5=60488505*h;
     var f1=95856640*h,f3=394240000*h,f4=206080000*h,
         f6=283115520*h;
     var g1=100061500*h,g3=375958000*h,g4=239527750*h,
         g5=18921870*h,g6=244823040*h;
     k1=fader(y0);
     for(j=1;n>j;j++){
         yn[j]=y0[j]+a1*k1[j];}
     k2=fader(yn);
     for(j=1;n>j;j++){
         yn[j]=y0[j]+b1*k1[j]+b2*k2[j];}
     k3=fader(yn);
     for(j=1;n>j;j++){
         yn[j]=y0[j]+c1*k1[j]+c2*k2[j]+c3*k3[j];}
     k4=fader(yn);
     for(j=1;n>j;j++){
         yn[j]=y0[j]+d1*k1[j]+d2*k2[j]+d3*k3[j]+d4*k4[j];}
     k5=fader(yn);
     for(j=1;n>j;j++){
         yn[j]=y0[j]+e1*k1[j]+e2*k2[j]+e3*k3[j]+e4*k4[j]+e5*k5[j];}
     k6=fader(yn);
     erromax=0.0;
     for(j=1;n>j;j++){
         delta=f1*k1[j]+f3*k3[j]+f4*k4[j]+f6*k6[j];
         yn[j]=y0[j]+delta;
         erromax=Math.max(erromax,Math.abs(g1*k1[j]+
               g3*k3[j]+g4*k4[j]+g5*k5[j]+g6*k6[j]-delta));}
     if(limite>erromax){
         for(j=0;n>j;j++) y0[j]=yn[j];
         return true;}
     else return false;
     };
/*
// -----
// Esempio di uso:
//
var prec=0.000001; // precisione errore.
var pasmed=20;     // passi necessari stimati.
// Supponendo che il tempo sia y0[3]...
// Debbono esistere anche i vettori k1...k6 e yn ed ym, tutti di 
// lunghezza uguale ad y0.
//
function moltipassi(){
    var j,k,hh,hhtp;
    hhtp=1/(pasmed*979292160.0);
    hh=hhtp;
    y0=[0,1.0,0,0];
    // Consento un numero sovrabbondante di passi... mille.
    for(j=1;1000>j;j++){
       if(cki(hh,prec)){
           y0[0]=j;
           for(k=0;4>k;k++) ym[k]=y0[k];
           // se avevo ridotto il passo per ragioni di precisione
           // lo riallungo lentamente, sperando di restare nei limiti
           // della precisione.
           hh=1.05*hh;
           // ma non supero mai il passo massimo ammesso.
           if(hh>hhtp) hh=hhtp;
           // Se sono vicino alla fine cerco di fare un passo opportuno
           // in modo da centrare esattamente l'istante voluto.
           if(y0[3]+hh*979292160.0 >2.000000001) hh=(2-y0[3])/979292160.0;
           }
       // Se supero la precisione ritento con passo dimezzato.
       else hh=0.5*hh;
       // Impongo la fine quando arrivo a due.
       if(y0[3]>=2.0) break;
       }
    return
    }
// 
//
*/