Capitolo 4 - la legge

Ora abbiamo quasi tutto quello che ci serve per costruire una simulazione della legge di gravitazione universale:

Il nostro obiettivo è simulare la traiettoria di un oggetto che si trova nel campo di attrazione gravitazionale di un altro oggetto virtuale che consideriamo inizialmente immobile e senza estensione, per semplicità. Per fare questo abbiamo bisogno di riscrivere le equazioni che determinano i valori istantanei del vettore accelerazione secondo la legge di Newton.

Per prima cosa abbiamo bisogno della distanza tra le due masse. Per trovarla, introduciamo due nuove variabili che indicano la posizione del centro di gravità.

float gx, gy; // coordinate del centro di gravità

Definiamo anche la costante di gravitazione universale che, nella nostra simulazione, potrà assumere valori arbitrari che comprenderanno anche la massa del corpo che genera il campo.

float G = 1000;

Per calcolare la distanza usiamo una funzione già disponibile tra quelle di base di Processing

float d = dist(x, y, gx, gy);

Come è facile intuire, i quattro parametri della funzione dist sono le coordinate dei due punti dei quali si vuole ricavare la distanza. E' utile prevedere una "istruzione di sicurezza" che, nel caso la distanza risulti esattamente nulla, eviti un blocco del programma dovuto ad una divisione per zero.

if (d==0) d = 1;

Notare la differenza tra il simbolo di uguaglianza logica "==" e quello dell'istruzione di assegnazione "=".

A questo punto siamo quasi pronti per ricavare le due componenti dell'accelerazione ricordando che è indipendente dalla massa del corpo in orbita, deve essere sempre diretta verso il centro di gravità e osservando che il suo vettore è dato dalla differenza vettoriale tra il vettore della posizione del centro di gravità e quello della posizione dell'oggetto orbitante. (v. fig. seguente, dove s è il vettore della posizione dell'oggetto, g è la posizione del centro di gravità, e a il vettore accelerazione)

In realtà in questo modo il vettore accelerazione avrà un modulo che è pari alla distanza tra oggetto e centro di gravità, per cui, per normalizzarlo e riportarlo ad un modulo unitario, dovremo dividerlo per la distanza stessa. Questo è il motivo dell'esponente 3 al quale viene elevata la distanza al denominatore della divisione:

ax = G * (gx - x) / pow(d, 3);

ay = G * (gy - y) / pow(d, 3);

La funzione pow, come è facile intuire, fornisce l'operazione di elevamento a potenza una volta dati base ed esponente.

Ed ecco il codice sorgente completo del nostro sistema di simulazione della legge di gravità universale:

// dichiarazione delle variabili

float x, y;       // due componenti per la posizione

float vx, vy;     // due componenti per la velocità

float ax, ay;     //accelerazione del punto

float gx, gy;     // coordinate del centro di gravità

float G = 1000;   // costante moltiplicativa

float t, dt;      // il tempo e l'incremento di tempo

float r = 20;     // raggio dell'oggetto

void setup()

{

  size(800, 600);

  t = 0;      // il contatore del tempo comincia da zero

  dt = 1;     // il delta-t viene inizializzato ad 1

  x = 100;    // x e y della posizione iniziale

  y = 200; 

  vx = 1;     // componenti x e y del vettore velocità

  vy = - 0.3;

  ax = 0.0;   // componenti iniziali dell'accelerazione

  ay = 0.1;

  gx = width / 2;     // coordinate del centro di gravità

  gy = height / 2;

}

void draw()

{

  //background(255);

  grafica();

  movimento();

  rimbalzo();

}

void grafica()

{

  fill(255, 0, 0);

  ellipse(x, y, 2*r, 2*r); // disegno l'ellisse

  textSize(15);                

  text("t: " + t, 10, 25);     // visualizzo il valore del tempo

}

void movimento()

{

  float d = dist(x, y, gx, gy);  // calcolo della distanza dal CdG

  if (d==0) d = 1;

  ax = G * (gx - x) / pow(d, 3); // 

  ay = G * (gy - y) / pow(d, 3);

  vx = vx + ax * dt;

  vy = vy + ay * dt;

  x = x + vx * dt;       // applico le equazioni del moto r.u.

  y = y + vy * dt;       // alla posizione

  t += dt;                  // incremento il tempo

}

void rimbalzo()

{

  if (x + r > width) // sto urtando la parete destra

  {

    vx = -vx;

    x = width - r;

  }

  if (x - r < 0)  // sto urtando la parete sinistra

  {

    vx = -vx;

    x = r;

  }

  if (y + r > height)  // sto urtando la parete inferiore

  {

    vy = -vy;

    y = height - r;

  }

  if (y - r < 0)  // sto urtando la parete superiore

  {

    vy = -vy;

    y = r;

  }

}

A partire da questo programma possiamo introdurre nuove caratteristiche del modello o simulare le traiettorie