//
//  Rungekutta ODE da integrare. 
//
opzionivalide=[0,1,2];
// _____________________________________________
//
//  In cosa consiste l'opzione 0
//
//  Problema del moto piano gravitazionale.
//
function variabili0(){
    return 5;
    }
function bloccodati0(){
    return 10;
    }
//
//  Il sistema di equazioni da integrare...
//
function sisder0(y,kr,ki){
    var c;
    y[kr] = 1;
    y[kr+1] = y[ki+3];
    y[kr+2] = y[ki+4];
    with(Math){
    c = - GperM*pow( y[ki+1]*y[ki+1] + 
                    y[ki+2]*y[ki+2] + 1.e-50, -1.5 );
          }
    y[kr+3] = c*y[ki+1];
    y[kr+4] = c*y[ki+2];
    }
//
//  Il vettore titolini[ ] contiene le didascalie mentre
//  il vettore valorini[ ] contiene i dati di inizializzazione
//
function promaschera0(){
    titolini[0]="Eccentricità dell'orbita";
    valorini[0]=0.01; 
    titolini[1]="Massa del sole [kg]";
    valorini[1]=1.9891e30;
    titolini[2]="Semiasse maggiore [m]";
    valorini[2]=1.50e11;
    titolini[3]="Anomalia eccentrica finale "+
        " ( ogni 3.14159 ossia PI, passa dal perielio all'afelio e viceversa "+
        " ossia effettua una semi_rivoluzione )";
    valorini[3]="2*PI";
    // alert("Esce da promaschera0()");
    }
//
//  Usa il vettore valorini[ ] per impostare i dati iniziali
//  del sistema di ODE.
//  Per effettuare l'inizializzazione si possono definire variabili
//  e function ausiliarie. In questo caso la variabile ecce e
//  il vettore analitico[ ] costruito tramite la function vero
//
var ecce=0.8;
var analitico=[0];
function vero(ecce,eta){
    // calcola la soluzione analitica per una data eccentricita'
    // ed una data anomalia eccentrica.
    var fattore,yv=[0];
    with(Math){
        fattore=sam*sqrt(sam/GperM)
        yv[0]=fattore*(eta-ecce*sin(eta));
        yv[1]=sam*(cos(eta)-ecce);
        yv[2]=sam*sin(eta)*sqrt(1-ecce*ecce);
        yv[3]=-sam*(sin(eta)/(1-ecce*cos(eta)))/fattore;
        yv[4]=sam*(cos(eta)*sqrt(1-ecce*ecce)/
                  (1-ecce*cos(eta)))/fattore;
        }
    return yv;
    }
//
var Gravuni=6.67428e-11; // costante di gravitazione universale.
var GperM=1;
var sam=1;   // semi asse maggiore
function sistemazione0(){
    var nobuono=true;
    if(valorini[0]>=0)if(1>valorini[0]){ 
            ecce=valorini[0];
            nobuono=false; };
    if(nobuono)alert("Errore l'eccentricità "+
       "deve stare tra 0 e 1. Usa il vecchio valore");
    GperM=valorini[1]*Gravuni;
    sam=valorini[2];
    analitico=vero(ecce,valorini[3]);
    tfin=analitico[0];
    // alert("Esce da sistemazione0 tfin="+(tfin/(24*3600))+" giorni");
    }
//
// Inizializza le derivate e specifica quali vuole stampare
// specificando gli elementi del vettore instampa[ ].
//
function inider0(y){
    var fattore=Math.sqrt(sam/GperM)
    y[0]= 0;
    y[1]= sam*(1-ecce); 
    y[2]= 0;
    y[3]= 0;
    y[4]= Math.sqrt((1+ecce)/(1-ecce))/fattore;
    instampa=[0,1,2,3,4]; // ossia stampa le variab. 0, 1,2,3 e 4.
    }
//
// Specifica il peso degli errori commessi. Ovviamente la scala 
// delle lunghezze è diversa da quella delle velocità.
//
function iniprec0(y){
    y[1]=1.e-9;
    y[2]=1.e-9;
    y[3]=1.e-3;
    y[4]=1.e-3;
    }
//
//  Calcoli finali di controllo. 
//  Vengono aggiunti dopo i dati dell'ultimo passo.
//
function confronti0(y){
    var j,np,ff="";
    if(0>nrestart) return ff;
    ff+=xlt+"br/>Soluzione esatta : ";
    for(j=0;instampa.length>j;j++) ff+=" y["+(instampa[j]%n_var)+"]= "+
           analitico[instampa[j]%n_var]+" ;";
    np=y[y.length-1];
    ff+=xlt+"br/>Utilizza il S.I. (  tempi in secondi, lunghezze in metri, "+
         "velocità in metri/secondo ) tfin ="+(y[np]/86400)+" [giorni].";
    return ff;
    }
//
// Per debuggare eventualmente i primi passi di calcolo.
//
function fapausa0(y){
     // stampa solo stringhe piu' lunghe di tre caratteri.
     // se il risultato è una stringa di 0 o 1 carattere
     // non verra' piu' chiamata. 
     if(totpassi>9) return "0";
     return "y[0]="+y[0]+xlt+"br/>";
     }
// _____________________________________________
//
//  In cosa consiste l'opzione 1
//
//  Problema del moto gravitazionale in tre dimensioni.
//
//  Costanti di servizio del problema 1:
//
var GperM1=1,Gravuni1=6.67428e-11,ecce1=0.02,analitico1=[0],sam1=1.5e11,aper1=0,incl1=0;
//
function variabili1(){
    return 7;
    }
function bloccodati1(){
    return 12;
    }
//
//  Il sistema di equazioni da integrare...
//
function sisder1(y,kr,ki){
    var c;
    y[kr] = 1;
    y[kr+1] = y[ki+4];
    y[kr+2] = y[ki+5];
    y[kr+3] = y[ki+6];
    with(Math){
    c = - GperM1*pow( y[ki+1]*y[ki+1] + 
                      y[ki+2]*y[ki+2] + 
                      y[ki+3]*y[ki+3] + 1.e-50, -1.5 );
          }
    y[kr+4] = c*y[ki+1];
    y[kr+5] = c*y[ki+2];
    y[kr+6] = c*y[ki+3];
    }
//
//  Il vettore titolini[ ] contiene le didascalie mentre
//  il vettore valorini[ ] contiene i dati di inizializzazione
//
function promaschera1(){
    titolini[0]="Eccentricità dell'orbita";
    valorini[0]=ecce1; 
    titolini[1]="Massa del sole [kg]";
    valorini[1]=1.9891e30;
    titolini[2]="Semiasse maggiore [m]";
    valorini[2]=sam1;
    titolini[3]="Angolo di rotazione dell'asse maggiore ( gradi e frazioni )";
    valorini[3]=aper1;
    titolini[4]="Inclinazione del piano orbitale ( gradi e frazioni )";
    valorini[4]=incl1;  
    titolini[5]="Anomalia eccentrica finale "+
        " ( ogni 3.14159 ossia PI, passa dal perielio all'afelio e viceversa "+
        " ossia effettua una semi_rivoluzione )";
    valorini[5]="PI";
    // alert("Esce da promaschera1()");
    }
//
//  Usa il vettore valorini[ ] per impostare i dati iniziali
//  del sistema di ODE.
//  Per effettuare l'inizializzazione si possono definire variabili
//  e function ausiliarie. In questo caso la variabile ecce1 e
//  il vettore analitico1[ ] costruito tramite la function vero1().
//
function vero1(ecce,eta,aper,incl){
    //
    // calcola la soluzione analitica per una data eccentricita'
    // ed una data anomalia eccentrica oltre che una data rotazione
    // dell'asse maggiore e una data inclinazione del
    // piano orbitale.
    //
    var xx,yy,fattore,yv=[0];
    with(Math){
        fattore=sam1*sqrt(sam1/GperM1)
        yv[0]=fattore*(eta-ecce*sin(eta));
        xx = [ sam1*(cos(eta)-ecce),
               sam1*sin(eta)*sqrt(1-ecce*ecce), 0.0];
        yy = xrota ( xx, aper, incl );
        yv[1] = yy[0];
        yv[2] = yy[1];
        yv[3] = yy[2];
        // alert("Posizioni "+yy);
        xx = [ -sam1*(sin(eta)/(1-ecce*cos(eta)))/fattore,
                sam1*(cos(eta)*sqrt(1-ecce*ecce)/
                  (1-ecce*cos(eta)))/fattore, 0.0 ];
        yy = xrota ( xx, aper, incl );
        yv[4] = yy[0];
        yv[5] = yy[1];
        yv[6] = yy[2];
        // alert("Velocita'"+yy);
        }
    return yv;
    }
//
function sistemazione1(){
    var nobuono=true;
    if(valorini[0]>=0)if(1>valorini[0]){ 
            ecce1=valorini[0];
            nobuono=false; };
    if(nobuono)alert("Errore l'eccentricità "+
       "deve stare tra 0 e 1. Usa il vecchio valore");
    GperM1=valorini[1]*Gravuni1;
    sam1=valorini[2]
    aper1=valorini[3];
    incl1=valorini[4];
    analitico1=vero1(ecce1,valorini[5],aper1,incl1);
    tfin=analitico1[0];
    // alert("Esce da sistemazione1 tfin="+(tfin/(24*3600))+" giorni");
    }
//
// Inizializza le derivate e specifica quali vuole stampare
// specificando gli elementi del vettore instampa[ ].
//
function inider1(y){
    var xx,yy,fattore=Math.sqrt(sam1/GperM1)
    y[0]= 0;
    xx=[ sam1*(1-ecce1), 0.0, 0.0 ];
    yy=xrota(xx,aper1,incl1);
    y[1]= yy[0]; 
    y[2]= yy[1];
    y[3]= yy[2];
    xx=[ 0.0, Math.sqrt((1+ecce1)/(1-ecce1))/fattore , 0.0 ];
    yy=xrota(xx,aper1,incl1);
    y[4]= yy[0];
    y[5]= yy[1];
    y[6]= yy[2];
    instampa=[0,1,2,3,4,5,6]; // ossia stampa tutte le variabili.
    }
//
function xrota(x,aper,incl){
    // Passa dalle coordinate del piano dell'orbita a quelle di calcolo.
    with(Math){
    var ca=cos(aper*PI/180);
    var sa=sin(aper*PI/180);
    var ci=cos(incl*PI/180);
    var si=sin(incl*PI/180);
        }
    return [ 
    (ci*(x[0]*ca - x[1]*sa)-x[2]*si),
    (x[0]*sa + x[1]*ca ),
    (si*(x[0]*ca - x[1]*sa)+x[2]*ci)];
    }
// 
//  Questa function non servirebbe ma e' utile averla a disposizione.
//
function arota(y,aper,incl){
    // inversa di xrota( ).
    var ca=Math.cos(aper*Math.PI/180);
    var sa=Math.sin(aper*Math.PI/180);
    var ci=Math.cos(incl*Math.PI/180);
    var si=Math.sin(incl*Math.PI/180);
    return [ (ca*(y[0]*ci + y[2]*si)+y[1]*sa) , 
             (y[1]*ca-sa*(y[0]*ci + y[2]*si)) , 
             (y[2]*ci-y[0]*si) ]; 
    }
//
// Specifica il peso degli errori commessi. Ovviamente la scala 
// delle lunghezze e' diversa da quella delle velocita'.
//
function iniprec1(y){
    y[1]=1.e-9;
    y[2]=1.e-9;
    y[3]=1.e-9;
    y[4]=1.e-3;
    y[5]=1.e-3;
    y[6]=1.e-3;
    }
//
//  Calcoli finali di controllo. 
//  Vengono aggiunti dopo i dati dell'ultimo passo.
//
function confronti1(y){
    var j,np,ff="";
    ff+=xlt+"br/>Soluzione esatta del problema tridimensionale : "+
        xlt+"br/>======= ";
    for(j=0;instampa.length>j;j++) ff+=" y["+(instampa[j]%n_var)+"]= "+
       analitico1[instampa[j]%n_var]+" ;";
    np=y[y.length-1];
    ff+=xlt+"br/>Utilizza il S.I. (  tempi in secondi, lunghezze in metri, "+
         "velocità in metri/secondo ) tfin ="+(y[np]/86400)+" [giorni].";
    return ff;
    }
//
// Per debuggare eventualmente i primi passi di calcolo.
//
function fapausa1(y){
     // stampa solo stringhe piu' lunghe di tre caratteri.
     // se il risultato è una stringa di 0 o 1 carattere
     // non verra' piu' chiamata. 
     return "0";
     }
// _____________________________________________
//
//  In cosa consiste l'opzione 2
//
//  Sistema binario in tre dimensioni.
//
//  Costanti di servizio del problema 2:
//
var ma2=1.9891e30,mb2=1.9891e27,Gravuni2=6.67428e-11,ecce2=0.01,sam2=7.6e11,analitico2=[0],aper2=0,incl2=2;
//
function variabili2(){
    return 13;
    }
function bloccodati2(){
    return 20;
    }
//
//  Il sistema di equazioni da integrare...
//
function sisder2(y,kr,ki){
    var dx,dy,dz,c,fat,rq,rr;
    y[kr] = 1;
    y[kr+1] = y[ki+7];
    y[kr+2] = y[ki+8];
    y[kr+3] = y[ki+9];
    y[kr+4] = y[ki+10];
    y[kr+5] = y[ki+11];
    y[kr+6] = y[ki+12];
    dx=y[ki+1]-y[ki+4];
    dy=y[ki+2]-y[ki+5];
    dz=y[ki+3]-y[ki+6];
    rq=dx*dx+dy*dy+dz*dz + 1.e-50;
    rr=Math.sqrt(rq);
    c= 1/(rq*rr); 
    fat= -c*mb2;
    y[kr+7] = fat*dx;
    y[kr+8] = fat*dy;
    y[kr+9] = fat*dz;
    fat= c*ma2;
    y[kr+10] = fat*dx;
    y[kr+11] = fat*dy;
    y[kr+12] = fat*dz;
    }
//
//  Il vettore titolini[ ] contiene le didascalie mentre
//  il vettore valorini[ ] contiene i dati di inizializzazione
//
function promaschera2(){
    titolini[0]="Eccentricità dell'orbita";
    valorini[0]=ecce2; 
    titolini[1]="Massa del primo corpo [kg]";
    valorini[1]=ma2;
    titolini[2]="Massa del secondo corpo [kg]";
    valorini[2]=mb2;
    titolini[3]="Semiasse maggiore [m]";
    valorini[3]=sam2;
    titolini[4]="Angolo di rotazione dell'asse maggiore ( gradi e frazioni )";
    valorini[4]=aper2;
    titolini[5]="Inclinazione del piano orbitale ( gradi e frazioni )";
    valorini[5]=incl2;  
    titolini[6]="Anomalia eccentrica finale "+
        " ( ogni 3.14159 ossia PI, passa dal perielio all'afelio e viceversa "+
        " ossia effettua una semi_rivoluzione )";
    valorini[6]="PI";
//  alert("Esce da promaschera2()");
    }
//
//  Usa il vettore valorini[ ] per impostare i dati iniziali
//  del sistema di ODE.
//  Per effettuare l'inizializzazione si possono definire variabili
//  e function ausiliarie. In questo caso la variabile ecce1 e
//  il vettore analitico1[ ] costruito tramite la function vero1().
//
function vero2(ecce,eta,aper,incl){
    //
    // calcola la soluzione analitica per una data eccentricita'
    // ed una data anomalia eccentrica oltre che una data rotazione
    // dell'asse maggiore e una data inclinazione del
    // piano orbitale.
    //
    var xx,yy,fattore,yv=[0];
    var fa=-mb2/(ma2+mb2);
    var fb=ma2/(ma2+mb2);
    //
    with(Math){
        fattore=sam2*sqrt(sam2/(ma2+mb2))
        yv[0]=fattore*(eta-ecce*sin(eta));
        xx = [ sam2*(cos(eta)-ecce),
               sam2*sin(eta)*sqrt(1-ecce*ecce), 0.0];
        yy = xrota2( xx, aper, incl );
        yv[1] = yy[0]*fa;
        yv[2] = yy[1]*fa;
        yv[3] = yy[2]*fa;
        yv[4] = yy[0]*fb;
        yv[5] = yy[1]*fb;
        yv[6] = yy[2]*fb;
        // alert("Posizioni "+yy);
        xx = [ -sam2*(sin(eta)/(1-ecce*cos(eta)))/fattore,
                sam2*(cos(eta)*sqrt(1-ecce*ecce)/
                  (1-ecce*cos(eta)))/fattore, 0.0 ];
        yy = xrota2( xx, aper, incl );
        yv[7] = yy[0]*fa;
        yv[8] = yy[1]*fa;
        yv[9] = yy[2]*fa;
        yv[10]= yy[0]*fb;
        yv[11]= yy[1]*fb;
        yv[12]= yy[2]*fb;
        // alert("Velocita' "+yy);
        }
    return yv;
    }
//
function sistemazione2(){
    var nobuono=true;
    if(valorini[0]>=0)if(1>valorini[0]){ 
            ecce2=valorini[0];
            nobuono=false; };
    if(nobuono)alert("Errore l'eccentricità "+
       "deve stare tra 0 e 1. Usa il vecchio valore");
    ma2=valorini[1]*Gravuni2;
    mb2=valorini[2]*Gravuni2;
    sam2=valorini[3];
    aper2=valorini[4];
    incl2=valorini[5];
    analitico2=vero2(ecce2,valorini[6],aper2,incl2);
    tfin=analitico2[0];
    // alert("Esce da sistemazione2 tfin="+(tfin/(24*3600))+" giorni");
    }
//
// Inizializza le derivate e specifica quali vuole stampare
// specificando gli elementi del vettore instampa[ ].
//
function inider2(y){
    var xx,yy,fa,fb,fattore=Math.sqrt(sam2/(ma2+mb2));
    y[0]= 0;
    xx=[ sam2*(1-ecce2), 0.0, 0.0 ];
    yy=xrota2(xx,aper2,incl2);
    fa=-mb2/(ma2+mb2);
    fb=ma2/(ma2+mb2);
    y[1]= fa*yy[0]; 
    y[2]= fa*yy[1];
    y[3]= fa*yy[2];
    y[4]= fb*yy[0];
    y[5]= fb*yy[1];
    y[6]= fb*yy[2];
    xx=[ 0.0, Math.sqrt((1+ecce2)/(1-ecce2))/fattore , 0.0 ];
    yy=xrota2(xx,aper2,incl2);
    y[7]= fa*yy[0];
    y[8]= fa*yy[1];
    y[9]= fa*yy[2];
    y[10]=fb*yy[0];
    y[11]=fb*yy[1];
    y[12]=fb*yy[2]
    instampa=[0,1,2,3,4,5,6,7,8,9,10,11,12]; // ossia stampa le var. del corpo b.
    }
//
function xrota2(x,aper,incl){
    // Passa dalle coordinate del piano dell'orbita a quelle di calcolo.
    with(Math){
    var ca=cos(aper*PI/180);
    var sa=sin(aper*PI/180);
    var ci=cos(incl*PI/180);
    var si=sin(incl*PI/180);
        }
    return [ 
    (ci*(x[0]*ca - x[1]*sa)-x[2]*si),
    (x[0]*sa + x[1]*ca ),
    (si*(x[0]*ca - x[1]*sa)+x[2]*ci)];
    }
// 
//  Questa function non servirebbe ma e' utile averla a disposizione.
//
function arota2(y,aper,incl){
    // inversa di xrota( ).
    var ca=Math.cos(aper*Math.PI/180);
    var sa=Math.sin(aper*Math.PI/180);
    var ci=Math.cos(incl*Math.PI/180);
    var si=Math.sin(incl*Math.PI/180);
    return [ (ca*(y[0]*ci + y[2]*si)+y[1]*sa) , 
             (y[1]*ca-sa*(y[0]*ci + y[2]*si)) , 
             (y[2]*ci-y[0]*si) ]; 
    }
//
// Specifica il peso degli errori commessi. Ovviamente la scala 
// delle lunghezze e' diversa da quella delle velocita'.
//
function iniprec2(y){
    y[1]=1.e-7;
    y[2]=1.e-7;
    y[3]=1.e-7;
    y[4]=1.e-7;
    y[5]=1.e-7;
    y[6]=1.e-7;
    y[7]=1.e-3;
    y[8]=1.e-3;
    y[9]=1.e-3;
    y[10]=1.e-3;
    y[11]=1.e-3;
    y[12]=1.e-3;
    }
//
//  Calcoli finali di controllo. 
//  Vengono aggiunti dopo i dati dell'ultimo passo.
//
function confronti2(y){
    var j,np,ff="";
    ff+=xlt+"br/>Soluzione esatta del problema tridimensionale : "+
        xlt+"br/>======= ";
    for(j=0;instampa.length>j;j++) ff+=" y["+(instampa[j]%n_var)+"]= "+
       analitico2[instampa[j]%n_var]+" ;";
    np=y[y.length-1];
    ff+=xlt+"br/>Utilizza il S.I. (  tempi in secondi, lunghezze in metri, "+
         "velocità in metri/secondo ) tfin ="+(y[np]/86400)+" [giorni] ="+
         (y[np]/(86400*365.25))+" [anni]";
    return ff;
    }
//
// Per debuggare, eventualmente, i primi passi di calcolo.
//
function fapausa2(y){
     // stampa solo stringhe piu' lunghe di tre caratteri.
     // se il risultato è una stringa di 0 o 1 carattere
     // non verra' piu' chiamata. 
     return "0";
     }
// _______________________________________
//