/* henon-heiles-sbab3c.c : file sorgente di un programma che consente di verificare numericamente che il metodo simplettico SBAB3C รจ del quarto ordine, controllando la conservazione dell'energia per il modello di Henon e Heiles. */ /* Piu' esattamente, questo file e' schema_henon-heiles-sbab3c.c, che e' VOLUTAMENTE INCOMPLETO rispetto al programma nel file henon-heiles-sbab3c.c, cosi' da facilitare la scrittura delle parti mancanti, come esercizio di programmazione. */ #include #include /* Definizione dei prototipi delle funzioni */ void corrector_motion(double x[], double y[], double omega[], double dt); void da_cartesiane_a_polari(double x, double y, double *rho, double *theta); void da_polari_a_cartesiane(double rho, double theta, double *x, double *y); double energia(double x[], double y[], double omega[]); void harm_osc_motion(double x[], double y[], double omega[], double dt); void perturb_motion(double x[], double y[], double dt); void sbab3(double x[], double y[], double omega[], double c_dt[], double d_dt[]); void sbab3c(double x[], double y[], double omega[], double c_dt[], double d_dt[], double coef_corr_dt); /* ### */ int main() { int i, numintervalli, ind_pot, pot10min, pot10max; double x[2], y[2], xiniz[2], yiniz[2], omega[2]; double energia, energiainiz; double c[2], d[2], coef_corr, Tf, dt, c_dt[2], d_dt[2], coef_corr_dt; char riga[200]; FILE *ifp; printf(" Leggo dal file di input i parametri e le condizioni iniziali\n\n"); /* Inizio delle operazioni di lettura dal file di input */ ifp = fopen("henon-heiles-sbab3c.inp" , "r"); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga le frequenze omega */ sscanf(riga, "%lf %lf", &omega[0], &omega[1]); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga di dati il valore dell'energia */ sscanf(riga, "%lf", &energiainiz); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga di dati i valori iniziali della seconda coppia di x e y */ sscanf(riga, "%lf %lf", &xiniz[1], &yiniz[1]); /* Pongo il valore iniziale della prima coordinata x uguale a 0 */ xiniz[0] = 0; /* Calcolo il primo momento y a partire dall'energia */ energia = energiainiz - omega[1] * ( pow(xiniz[1],2) + pow(yiniz[1],2) ) / 2 + pow(xiniz[1],3) / 3; yiniz[0] = sqrt( 2 * energia / omega[0] ); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga il tempo finale di integrazione */ sscanf(riga, "%lf", &Tf); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga gli esponenti (in base 10) del numero di intervalllini */ sscanf(riga, "%d %d", &pot10min, &pot10max); /* Fine delle operazioni di lettura dal file di input */ fclose(ifp); /* Pongo i valori degli elementi dei vettori c e d uguali alle costanti relative al metodo SBAB3 (vedasi articolo di Laskar e Robutel del 2001 su Celestial Mechanics and Dynamical Astronomy) */ c[0] = 0.5 - sqrt(5.) / 10.; c[1] = sqrt(5.) / 5.; d[0] = 1. / 12.; d[1] = 5. / 12.; /* Pongo il valore del coefficiente di correzione uguale all'opportuna costante del metodo SBAB3C (vedasi il suddetto articolo di Laskar e Robutel del 2001) */ coef_corr = ( 13. - 5. * sqrt(5) ) / 288.; /* Stampo sullo schermo alcune righe che spiegano come interpretare i risultati che verranno visualizzati in seguito */ printf(" Numintervalli Delta E\n"); printf("-----------------------------\n"); for(ind_pot=pot10min; ind_pot<=pot10max; ind_pot++) { /* Pongo i valori delle coordinate canoniche uguali ai corrispondenti valori iniziali letti dal file di input */ for(i=0; i<2; i++) { x[i] = xiniz[i]; y[i] = yiniz[i]; } /* Pongo numintervalli uguale a 10 elevato alla i */ numintervalli = pow(10 , ind_pot); /* Pongo l'intervallino di integrazione uguale al rapporto tra Tf e numinitervalli */ dt = Tf / numintervalli; /* Calcolo i valori degli elementi dei vettori c_dt e d_dt in modo che corrispondano ai valori dei passi di integrazione che vengono effettuati (lungo il flusso delle due parti in cui viene divisa la Hamiltoniana) durante l'applicazione del metodo SBAB3 */ for(i=0; i<2; i++) { c_dt[i] = c[i] * dt; d_dt[i] = d[i] * dt; } /* Calcolo il valore di coef_corr_dt in modo che sia uguale al valore dei passi di integrazione che vengono effettuati lungo il flusso del termine di correzione, durante l'applicazione del metodo SBAB3C */ coef_corr_dt = - coef_corr * 0.5 * dt * dt * dt; /* Calcolo i valori delle coordinate canoniche al tempo Tf applicando numintervalli volte il metodo SBAB3C */ for(i=0; i