Calcolo della distanza geodetica tra due punti della superficie terrestre

Diversi anni fa, quando lavoravo al CNR, insieme a Bruno Fornari, un caro collega ed amico purtroppo scomparso, scrivemmo un programmino di "servizio" che utilizzammo per calcolare matrici di distanze terrestri (geodetiche). Il programma è veramente semplice, più interessante è capire come si calcola la distanza tra due punti di una superficie sferica. Devo riconoscere che, al momento, ho avuto la sensazione che i miei studi scientifici non fossero serviti a molto, perché non avevo la più pallida idea di come si facesse! (Gli studi dell'OCSE-Pisa dicono il vero...) Per fortuna, i consigli di Bruno, fisico e scienziato di prim'ordine, mi hanno fornito la formula giusta.
Per il nostro scopo, non ci serviva una grande precisione, per cui l'assunto iniziale era che la terra fosse una sfera perfetta (in realtà non lo è, vi rimando ai link sotto per gli approfondimenti). Dunque, osservando la figura in alto, diciamo che, in base alla trigonometria sferica (teorema di Eulero), tra i lati a, b e p del triangolo sferico ABP vale la relazione:

cos p = cos a cos b + sen a sen b cos φ

Ora, dette lat(A), lon(A), lat(B), lon(B), la latitudine e la longitudine dei punti A e B e, considerato che:

• a = 90° - lat(B)
• b = 90° - lat(A)
• φ = lon(A) - lon(B)

abbiamo tutti i dati per calcolare la lunghezza del lato p considerando il raggio della Terra approssimabile a R = 6371 km.

Ecco dunque il codice in linguaggio ANSI C:

dis_geod.h

#include <math.h>
/* Questa funzione calcola la distanza tra due punti 
sulla superficie terrestre, date le coordinate in
latitudine e longitudine espresse in
gradi decimali */
double disgeod (double latA, double lonA,
                double latB, double lonB)
{
      /* Definisce le costanti e le variabili */
      const double R = 6371;
      const double pigreco = 3.1415927;
      double lat_alfa, lat_beta;
      double lon_alfa, lon_beta;
      double fi;
      double p, d;
      /* Converte i gradi in radianti */
      lat_alfa = pigreco * latA / 180;
      lat_beta = pigreco * latB / 180;
      lon_alfa = pigreco * lonA / 180;
      lon_beta = pigreco * lonB / 180;
      /* Calcola l'angolo compreso fi */
      fi = fabs(lon_alfa - lon_beta);
      /* Calcola il terzo lato del triangolo sferico */
	  p = acos(sin(lat_beta) * sin(lat_alfa) + 
        cos(lat_beta) * cos(lat_alfa) * cos(fi));
      /* Calcola la distanza sulla superficie 
      terrestre R = ~6371 km */
      d = p * R;
      return(d);
}

Ed ecco il codice per il programma di interfaccia utente:

prodisge.c

#include <stdio.h>
#include "dis_geod.h"
double disgeod(double latA, double lonA, 
              double latB, double lonB);
void main (void)
{
float latA, lonA, latB, lonB, distanza;
printf("\n Inserisci la latitudine del punto A :\t");
scanf ("%f", &latA);
printf("\n Inserisci la longitudine del punto A :\t");
scanf ("%f", &lonA);
printf("\n Inserisci la latitudine del punto B :\t");
scanf ("%f", &latB);
printf("\n Inserisci la longitudine del punto B :\t");
scanf ("%f", &lonB);
distanza =  disgeod(latA, lonA, latB, lonB);
printf("\n La distanza fra A e B e' : %f  km\n", distanza);
}

Ed ecco infine il risultato del calcolo della distanza fra Roma (lat: +41.91;, lon: +12.45) e Milano (lat: +45.48 lon: +09.18)

Se provate il calcolo con altri programmi, potrete avere risultati un poco differenti, questo è dovuto alle approssimazioni usate per il raggio della Terra, per il valore di Π e per la conversione dei gradi sessagesimali in decimali.
Ovviamente, al programma originale fu aggiunto il codice necessario per costruire le matrici di distanze, che qui, per brevità, ho tagliato.

Download: (compilato per win32) prodisge.exe

Conclusioni:

Ormai, che i GPS ci parlano amichevolmente in auto, e Google Earth ci da tutte le informazioni che vogliamo con pochi click, programmini del genere sembrano della preistoria. Ma, a furia di avere sempre la "pappa pronta", non ci annichiliremo un po' troppo le meningi? Trovo che ogni tanto sarebbe bene fare qualche ripasso, tanto per capire anche come funziona ciò che si utilizza passivamente!

Riferimenti ed approfondimenti:

• Vittorio Villasmunta: Corso Basico di meteorologia
• Mauro Bertolini: Mini corso di trigonometria sferica
• La geometria sulla sfera: Un approccio elementare alla geometria non euclidea
• Convert Degrees Minutes Seconds to Decimal Degrees
• Federico Peiretti: Matematica e... la misura del mondo
• Donata Allegri: La Terra si sta allargando
• Navigazione Ortodromica e lossodromica
• Ed Williams: Great Circle Calculator
• Lega Navale Italiana - Richiami di geografia e sulle coordinate geografiche