La formula di Bowring

Il primo modello matematico di riferimento adatto a descrivere la superficie terrestre conrisultati apprezzabili che supera l’ingenua idea di attribuire alla terra una forma perfettamentesferica e che non è particolarmente complesso nella sua espressione matematica è statoindividuato fin dalla seconda metà del XVII secolo nell’ellissoide di rotazione.

Lo scopo di questo lavoro è quello di enunciare, dimostrare ed applicare la formula diBowring che è spesso utilizzata nei sistemi che implementano il modello dell’ellissoide dirotazione.

1.1 Il modello matematico : definizioni

Consideriamo un sistema cartesiano ortonormaleOxyzi cui versori identificativii,j,kformano una terna destrorsa.

Supponiamo di tracciare nel pianoOxzuna ellisse centrata inO con semiasse maggiore dilunghezzaa > 0 rivolto verso l’assex e semiasse minore di lunghezzab > 0 rivolto all’assez.

La superficie matematicaΣa,b che si ottiene ruotando questa ellisse attorno all’assezsi diceellissoide di rotazionee la sua equazione cartesiana è

x2 + y2

a2 + z2

b2 − 1 = 0 (1.1.1)

L’ellisse di semiassia,b contenuta nel pianoOxzsi diceellisse meridiana, mentre l’assezattorno cui ruota l’ellisse meridiana si diceasse polare.

L’ equatore Eè il luogo geometrico dei punti appartenenti al pianoOxy, dettopianoequatoriale, che hanno la medesima distanzaa dal centroO

E := x,y, 0 ∈ Oxy : x2 + y2 = a2 (1.1.2)

Con questa terminologia il semiasse maggiorea prende il nome disemiasse equatoriale, mentreil semiasse minoreb si dicesemiasse polare.

Si dicemeridianouna qualsiasi sezione piana (ellisse) ottenuta secando l’ellissoide con unpiano contenente l’asse polare.

Si diceparallelouna qualsiasi sezione piana (circonferenza) ottenuta secando l’ellissoidecon un piano parallelo al piano equatoriale.

Il punto di intersezioneN tra la superficie dell’ellissoide di rotazione e l’assezdiretto nelverso positivo si dicepolo nord, mentre l’intersezioneS ottenuta nella direzione opposta si dicepolo sud.

Le loro coordinate valgono

N = 0,0,b

S = 0,0,−b

(1.1.3 - I)

(1.1.3 - II)

Si dicesemimeridianopassante per un puntoP ∈ Σa,b la linea ottenuta come intersezione dellasuperficie ellissoidica con il semipiano formato dall’asse polare e passante per il puntoPassegnato; in alternativa si può definirla in modo del tuttoequivalente come la semiellissecompresa tra i poli nord e sud e passante per il punto datoP.

E’ invalso l’uso di chiamare il semimeridiano più semplicemente soltanto come ‘meridiano’,mentre il semimeridiano passante per l’antipodo del generico puntoP si diceantimeridiano.

Adottando questa convenzione definiamo inoltre ilmeridiano fondamentalela semiellisse

passante per i poli nord e sud e per il punto di coordinate1,0,0.L’ellissoide risulta essere completamente determinato quando sono noti i valoria,b dei

semiassi, oppure quando è noto soltanto uno dei due semiassie uno dei seguenti parametriadimensionali

schiacciamento : α = a − ba

prima eccentricità : e = a2 − b2

a2

seconda eccentricità : e′ = a2 − b2

b2

(1.1.4 - I)

(1.1.4 - II)

(1.1.4 - III)

1.2 Sistemi di coordinate globali

Nel modello appena descritto è naturale adottare tre tipologie differenti di sistemi diriferimento che hanno l’origine in comune proprio nel baricentro geometrico dell’ellissoide(puntoO).

Si badi comunque che per ora questi sistemi non forniscono nèrappresentano una definizioneideale di sistema di riferimentosolidalecon la terra sebbene ne costituiscano una importanteanticipazione, tuttavia nel prossimo paragrafo risolveremo la questione spiegando in modoestremamente sintetico quali sono stati i criteri comunemente scelti su scala globale per giungeread una definizione universale.

1. Coordinate cartesiane ellissocentricheE’ proprio il sistema cartesianoOxyzintrodotto nella sezione precedente e definisce in

maniera univoca una posizione plano-altimetrica in ambiente 3D. Per comodità ripetiamo le suecaratteristiche fondamentali: l’origineO coincide con il centro dell’ellissoide, l’asseOx è rivoltoverso il meridiano fondamentale e giace sul piano equatoriale, l’asseOzè rivolto verso l’assepolare (cioè di rotazione dell’ellisse meridiana) e l’asseOy è tale che la ternaOxyzsia di tipodestrorsa.

2. Coordinate angolari ellissocentricheE’ una coppia di angoliψ,ω che definisce in modo univoco un puntoP ∈∑a,b

dellasuperficie ellissoidica.

Latitudineψ = angolo compreso tra il raggio vettoreOP e il piano equatoriale, contato versonord (latitudine N o positiva) o verso sud (latitudine S o negativa) e limitato agli intervalli

0° ≤ ψ ≤ 90°N

0° ≤ ψ ≤ 90°S

Longitudineω = angolo compreso tra il piano del meridiano passante per il punto P ∈∑a,be il

piano del meridiano fondamentale , contato verso est (longitudine E o negativa) o verso ovest(longitudine W o positiva) e limitato agli intervalli

0° ≤ ω ≤ 180°E

0° ≤ ω ≤ 180°O

3. Coordinate angolari ellissoidiche

Dal punto di vista planimetrico si tratta di una coppia di angoli ϕ,ω atta a definire in modounivoco un puntoP ∈∑a,b

della superficie ellissoidica.Latitudineϕ = angolo compreso tra la normale ellissoidica perP e il piano equatoriale,

contato verso nord (latitudine N o positiva) o verso sud (latitudine S o negativa) e limitato agliintervalli

0° ≤ ϕ ≤ 90°N

0° ≤ ϕ ≤ 90°S

L’angoloω è la longitudine diP e coincide perfettamente con la definizione data al caso delsistema di coordinate geografico ellissocentrico.

Per effettuare una osservazione tridimensionale completadella posizione di un punto , sidefinisce anche laquota ellissoidica hche è la distanza di un puntoPh dalla sua proiezionePsull’ellissoide misurata lungo la normale ellissoidica.

1.3 - Sistemi solidali alla terra

La forma reale della terra è molto complessa perchè scendendo nel dettaglio si presentainfatti molto frastagliata sulla terraferma essendo costituita da rilievi montuosi, zone pianeggiantie collinari, inoltre essa è anche in continua evoluzione dinamica nel tempo perchè ad esempio lemasse oceaniche sono sottoposte a deformazioni giornaliere (maree) per effetto dell’attrazionegravitazionale combinata della luna e del sole.

La rappresentazione più ‘naturale’ della superficie terrestre nasce idealmente come ilprolungamento totale della superficie marina del pianeta al netto delle variazioni mareali e vienedettageoide.

Essa viene utilizzata correntemente per la rilevazione delle quote altimetriche in numeroseapplicazioni topografiche, ma la sua equazione rimane nonostante ciò ancora troppo complicata,da qui la necessità di approssimare la superficie terrestreassimilandola per comodità a quella diun ellissoide di rotazione (v. fig. sottostante).

Fig. (1.3.1)

Nel paragrafo 1.1 abbiamo fornito le prime definizioni a carattere matematico dell’ellissoide,tuttavia non abbiamo ancora specificato quantitativamente la forma dell’ellissoide (semiassemaggiore e appiattimento) e il suo orientamento spaziale.

Per quanto riguarda l’orientamento nella geodesia satellitare (sistemi di tipo GPS eGLONASS) attualmente si utilizzano sistemi di misurazionecartesianiOxyzin cui

l’origine O viene fatta coincidere con il centro di massa del pianeta , l’asse polare zcoincide con l’asse di rotazione terrestre e gli assi x,y completano la terna destrogiratrovandosi collocati sul piano equatoriale, con l’asse x rivolto verso il meridiano fondamentalepassante per l’osservatorio di Greenwich.

Sistemi di questo tipo vengono detti sistemiECEF (Earth-Centered-Earth-Fixed) in quantola terna cartesiana è perfettamente solidale alla terra, comunque va detto che la rilevazione delcentro di massa terrestre e la collocazione dell’asse di rotazione terrestre sono stabiliti in base acriteri puramente convenzionali, sebbene siano supportati da misurazioni satellitari di altissimaprecisione.

Il sistema ECEF maggiormente utilizzato in questi anni che èstato implementato nella retesatellitare GPS è ilWGS84 (World Geodetic System ’84) in cui i parametri dell’ellissoide sonostati definiti nel 1984 in questo modo:

a = 6.378.137 m.

α = 1298,2572221

Gli altri tre sistemi più importanti ancora utilizzati al giorno d’oggi sono l’Hayford (1924 -Associazione internazionale di geodesia) con parametri

a = 6.378.388 m.

α = 1297

il Bessel (calcolato nel 1841) di parametri

a = 6.377.397,155m.

α = 1299,152

e l’ellissoide internazionale (1967) i cui parametri sono stati riscontrati grazie alle osservazionieffettuate dai satelliti artificiali

a = 6.378.160 m.

α = 1298,017

ed è quello che approssima meglio di tutti il geoide.In base alle specifiche ECEF non converrà più parlare di sistemi ellissocentrici anche se da

un punto di vista geometrico la denominazione continua ad essere formalmente corretta, maparleremo invece disistemi geocentrici.

In tal caso le denominazioni utilizzate nella lessicografia topografica diventano :

Coordinate cartesiane geocentriche x,y,z Coordinate cartesiane ellissocentriche

Coordinate geografiche ϕ,ω Coordinate angolari ellissoidicheIl parametroϕ viene detto latitudine geografica (o ellissoidica) di un punto collocato

sull’ellissoide di riferimento.

Coordinate angolari geocentriche ψ,ω Coordinate angolari ellissocentricheIl parametroψ viene detto latitudine geocentrica di un punto collocato sull’ellissoide di

riferimento.

1.4 - Trasformazioni di coordinate : teoremi

Theorem 1.4.1 - Trasformazione diretta da coordinate geografiche a cartesiane geocentriche

ϕ,ω,h x,y,z

x = N + hcosϕcosω

y = N + hcosϕsinω

z = N1 − e2 + hsinϕ

(1.4.1 -I)

(1.4.1 -II )

(1.4.1 -III )

dove

N = a 1 − e2 sin2ϕ − 12 (1.4.1 -IV)

Theorem 1.4.2 - Trasformazione inversa da coordinate cartesiane geocentriche a geografiche

x,y,z ϕ,ω,h

ϕ = arctanq

ω = arctanyx ∗

h =x2 + y2

cosϕ − N

(1.4.2 -I)

1.4.2 -II

(1.4.2 -III )

dove q è soluzione di x2 + y2 q − ae2q

1+ 1−e2 q2− z = 0

Nota (*) : l’arcotangente ha codominio in− π2 ,+ π

2sex = 0 , y = 0 ω = 0° (per convenzione)sex = 0 : y ≷ 0 ω = ± 90°sex > 0 ω = arctany

x

sex < 0 , y ≥ 0 ω = π + arctanyx

sex < 0 , y < 0 ω = −π + arctanyx

Theorem 1.4.3 - Trasformazione diretta da coordinate angolari geocentriche a cartesianegeocentriche. E’ una parametrizzazione della superficie ellissoidica.

ψ,ω x,y,z

x = Rcosψcosω

y = Rcosψsinω

z = Rsinψ

(1.4.3 -I)

(1.4.3 -II )

(1.4.3 -III )

dove

R = a 1 − e2

1 − e2 cos2ψ (1.4.3 -IV)

Theorem 1.4.4 - Relazione matematica tra la latitudine geograficaϕ e la latitudinegeocentricaψ valida solo per i punti dell’ellissoide

tanψ = 1 − e2 tanϕ (1.4.4)

Theorem 1.4.5 - Formula della latitudine geografica di Bowring

ϕ = arctan z+ e′2bsin3θs− e2acos3θ

dove

s = x2 + y2

θ = arctan zasb

e′2 = a2 − b2

b2 (seconda eccentricità)

1.5 - Deduzione delle formule di trasformazione

In questo paragrafo deduciamo tutte le formule sulle trasformazioni da un sistema dicoordinate all’altro riunendo in un unico schema sequenziale tutti i vari passaggi logicifondamentali.

Consideriamo dunque una sezione di ellissoide ottenuta perintersezione con il pianopassante perP = Pr,z contenente l’asse polareOz, come indicato in figura.

Fig. 1.5.1

La direzioner individua l’asse contenuto nel piano equatoriale che ha origine inO e puntaverso la proiezionePr di P sul medesimo piano.

Le due coordinater,zdi P non sono libere perchè sono legate attraverso l’equazione delmeridiano passante perP

r2

a2 + z2

b2 − 1 = 0 (1.5.1)

che identifica come già detto una ellisse di semiassia,b.La normale ellissoidica (non normalizzata) al puntoP nel pianoOrzsi ottiene calcolando il

gradiente della equazione (1.5.1)

n = ra2 , z

b2 (1.5.2)

Il punto P ′ è l’intersezione tra la retta normale alla superficie passante inP e l’asseOr. Perdeterminarne esplicitamente le coordinate dobbiamo quindi parametrizzare la retta normale conuna nuova variabile (che chiamiamo ad esempiot) e annullare la seconda componente (z) :

z+ t∗ zb2 = 0 t∗ = −b2

OP′= r + t∗ r

a2 = r 1 − b2

a2 = re2 (1.5.3)

Dalla relazione geometrica

OP′+ P′Pr = OPr (1.5.4)

si ricava immediatamente

P′Pr = r − re2 = r1 − e2 (1.5.5)

Poichè l’altezzazdi P è funzione della sua ascissar secondo la relazione (1.5.1) , alloraricaviamo facilmente anche

z = ba a2 − r2 = 1 − e2a2 − r2 (1.5.6)

e grazie al teorema di Pitagora applicato al triangolo rettangoloΔP′PrP possiamo finalmentededurre la lunghezzad dell’ipotenusaP′P sommando i quadrati dei cateti ed estraendo la radice

d = r21 − e22 + 1 − e2a2 − r2 =

= 1 − e2a2 − r2e2 (1.5.7)

Tale relazione che esplicita il legame esistente tra i parametri r,d può avere una sua importanzaautonoma, perciò in certi casi è bene servirsi anche della sua inversa

r = 1e a2 − d2

1 − e2 (1.5.8)

Introduciamo ora l’angoloϕ = PrP′P che è la latitudine geografica diP e sfruttiamo di nuovo la(1.5.4) ottenendo così

dcosϕ = r1 − e2 r =dcosϕ1 − e2 (1.5.9)

Sostituendo lar di (1.5.9) in (1.5.7) si arriva ad avere una equazione nell’incognitad2

d2 = 1 − e2 a2 −d2 cos2ϕ

1 − e22 e2

che risolta estraendone la radice quadrata dà

d =a1 − e2

1 − e2 sin2ϕ= N1 − e2 (1.5.10)

doveN è il parametro indicato in (1.4.1 - IV).Per ottenere la relazione (1.4.1 - III) dell’altezzaz in funzione diϕ (a quota ellissoidica

nulla) è sufficiente eseguire il prodottodsinϕ , mentre se vogliamo ricavare le restanti equazioniparametriche inx,y dobbiamo prima riscrivere nuovamente la (1.5.9) per mezzo della relazione(1.5.10)

r =N1 − e2cosϕ

1 − e2 = Ncosϕ (1.5.11)

e infine applicare ovviamente le equazioni

x = r cosω

y = r sinω

in modo tale che anche le (1.4.1 - I, II) restano identicamente soddisfatte.Le corrispondenti formule planoaltimetriche si ottengonoevidentemente aggiungendo ad il

valoreh della corrispondente quota ellissoidica, per cui il teorema (1.4.1) resta così dimostrato.Supponiamo ora di invertire le formule del teorema (1.4.1) dando per scontato per un

momento di conoscere il valore esatto della latitudineϕ.Sfruttando ora le prime due equazioni del teorema (1.4.1) è molto facile ricavare il valore di

h in funzione delle sole variabili cartesiane equatorialix,y e diϕ, infatti banalmente si trova

h =x2 + y2

cosϕ − N

che è proprio la terza equazione del teorema (1.4.2).Dividendo poi la componentey per la componentex si perviene infine anche alla seconda

relazione di (1.4.2)

tanω =yx , x ≠ 0

x = 0 : y ≷ 0 ω = ± 90°

a cui però è necessario prestare attenzione all’uso della funzione arcotangente per determinare lagiusta collocazione diω in uno solo dei quattro quadranti del piano equatorialeOxy(v. nota nellasez. 1.4 , pag. 6).

Molto più problematica da ottenere è invece la soluzione esplicita della latitudine geograficaϕ in funzione delle sole componenti cartesiane. In letteratura sono riportati numerosiprocedimenti matematici per ottenereϕ, alcuni dei quali sono algoritmi iterativi come quello checi accingiamo a mostrare.

Tuttavia una soluzione in formula chiusa esiste ed è stata scoperta da Bowring (v. Formula diBowring, theorem 1.4.5) e dimostreremo anch’essa alle pag.12 - 13.

Partiamo allora dalla terza equazione di (1.4.1) e sostituiamo il parametroh con la formula(1.4.2 - III) che è stata appena riscoperta poco più sopra, poi trasformiamo il tutto eliminando lefunzioni seno e coseno favorendo la comparsa della tangentedi ϕ :

z =x2 + y2

cosϕ − ae2

1 − e2 sin2ϕ12

sinϕ =

= x2 + y2 tanϕ −ae2 tanϕ

1 + 1 − e2 tan2ϕ (1.5.12)

ed effettuando la sostituzioneq = tanϕ ritroviamo l’equazione (1.4.2 - I) che incidentalmente èanche una formula alternativa per calcolare l’altezza cartesianaz in funzione delle coordinatex,ye della latitudineϕ, quindi dopo aver trovato la sua soluzioneq∗ si calcolaϕ attraversol’arcotangente.

Per impostare correttamente l’algoritmo iterativo occorre pertanto studiare preliminarmentel’andamento della funzione

zq = qs− ae2q

1 + 1 − e2q2, con s= x2 + y2 (1.5.12 bis)

sapendo che i parametria,e,ssono positivi e che la radice quadrata a denominatore è ovunque

definita in quanto

1 − e2 = 1 − a2 − b2

a2 = b2

a2 > 0

Dunque la funzionezq è sempre continua inR, inoltre è dispari (e quindi passantenell’origine).

La derivata prima e seconda dizq sono

z′q = s− ae2

1 + 1 − e2q232

z′′q =3ae21 − e2q

1 + 1 − e2q252

La derivata prima è positiva soltanto quando vale la condizione

q2 > 11 − e2

ae2

s

23− 1 (1.5.13)

mentre la derivata seconda indica chiaramente che la funzione volge la concavità verso l’altosolo perq ≥ 0.

La funzionezq è sempre crescente quandoae2 < s, altrimenti in caso contrario quandoae2 > ssi ha un intero intervallo simmetrico centrato nell’origine in cuizq ha andamentodecrescente.

In figura (1.5.14) sono riportati tre grafici relativi allafunzionezq nei casiae2 ⋚ ssottol’ipotesi che i parametrie2 = 0.3 ,s = 6 siano mantenuti costanti.

Fig. 1.5.14

Come si vede nel riquadro centrale quando vi è perfetta uguaglianzas = ae2 , allorazqpresenta un flesso a tangente orizzontale inq = 0.

Un eccellente algoritmo iterativo che si presta particolarmente bene a trovare la radiceq∗

della equazionez− zq = 0 è il notometodo di Newtonche ha, com’è noto dalla teoria, ordinedi convergenza pari a due e formalmente si scrive

qk+1 = qk −zqk − z

z′qk

Applicandolo specificatamente al caso della funzione (1.5.12) il metodo si esplicita nel modoseguente

A = 1 + 1 − e2qk2

qk+1 = qk − A qks− z A − qkae2

s A3 − ae2 (1.5.15)

Per quanto riguarda l’innesco del metodo, dobbiamo distinguere due casi particolari1. ae2 ≤ s : è il caso più semplice, infatti basta cercare un valore iniziale q0 tale che sez > 0

allorazq0 > ze se invecez < 0 allorazq0 < z. Naturalmente sez = 0 la soluzione èimmediatamenteq = 0, cioèϕ = 0.

2. ae2 > s: la funzionezq non è globalmente monotòna e fissando un valorezpotrebbeesserci purtoppo più di una controimmagine (fino ad un massimo di tre), comeeffettivamente si evince dal terzo riquadro della fig. (1.5.14). In tal caso il criterio di scelta èpuramente convenzionale: per esempio possiamo stabilire di scartare le soluzioni che hannosegno opposto alla quotaz. In altri termini, sez ≷ 0 alloraq∗ ≷ 0, sicchè il valore diinnescoq0 può essere cercato come indicato al punto 1.

Il caso 2 ha una interessante interpretazione geometrica, infatti attesta l’esistenza di uninsieme di punti appartenenti al piano meridiano che godonodella proprietà di trovarsiall’incrocio delle direzioni di più normali ellissoidiche. Nella figura sottostante è riportato unesempio in cui si vede che il puntoP è posizionato lungo tre direzioni normali, di cui unasoltanto è positiva mentre le altre due sono negative.

Fig. (1.5.16)

Il criterio di arresto del metodo consiste nel confrontare il valore assoluto della differenza tradue approssimazioni consecutive diq con un prefissato valore positivo sufficientemente piccolo(la tolleranza), in modo tale che l’algoritmo termini non appena questa differenza diventa

inferiore.Il grado di precisione del metodo è quindi legato sostanzialmente all’ordine di grandezza

della tolleranza e deve essere scelto in conformità al tipo di applicazione che si desideraeffettuare, ad esempio nel sistema WGS84 commettere un errore massimo posizionale di 1 mt.lungo un meridiano significa accettare una tolleranza massima angolare pari a

arctan 16.378.137

≅ 1.567e− 007rad

cioè in definitiva pari a circa 0.032′′ (secondi d’arco).Passiamo ora alle coordinate geocentriche e cerchiamo di studiare nel dettaglio le relazioni

matematiche esistenti con i sistemi di coordinate precedenti.Intanto ricordando le definizioni del paragrafo 1.2 e facendo riferimento alla figura (1.3.1)

che fornisce una visualizzazione astratta dei sistemi geografici globali, sappiamo che lalongitudineω è assolutamente identica sia al caso del sistema geocentrico che al caso di quelloellissoidico.

Indichiamo conR la distanza di un puntoPr,z appartenente alla superficie ellissoidicadall’origine O del sistema geocentrico, allora è facile calcolare la lunghezza delle proiezioniz,rdi P rispettivamente sull’asse polare e sul piano equatoriale

z = |OPz| = Rsinψ

r = |OPr | = Rcosψ

(1.5.17 - I)

(1.5.17 - II)

Adesso andiamo a sostituirle nell’equazione del meridiano(1.5.1) ricavando così la lunghezza(sempre positiva) del raggio vettoreOP in funzione della latitudine geocentrica

R = a 1 − e2

1 − e2 cos2ψ (1.5.18)

viceversa, se disponiamo del valore della distanzaR, possiamo invertire la formula e ricavare lalatitudine (positiva e negativa)

ψ = ±arctan pe2 − p

, dovep = 1 − R2

a2 (1.5.19)

però dobbiamo ovviamente prestare attenzione al caso in cuiR = b = a 1 − e2 poichè ildenominatore dentro la radice si annulla e gli angoli valgono naturalmenteψ = ±90°.

Conoscendo la lunghezzaR del raggio vettore si deducono facilmente le formuleparametriche dell’ellissoide (1.4.3 - I, II, III) nelle variabili ψ,ω.

Confrontando poi la (1.5.11) con la (1.5.17 - II) si ottiene

Rcosψ = Ncosϕ (1.5.20)

e confrontando la (1.4.3 - III) con la (1.4.1 - III) a quota ellissoidica nulla si ha invece

Rsinψ = N1 − e2sinϕ (1.5.21)

per cui dividendo membro a membro la (1.5.21) con la (1.5.20)ricaviamo finalmente larelazione del teorema (1.4.4).

Adesso passiamo invece a dimostrare l’importante formula di Bowring (1.4.5) che fornisceesplicitamente il valore della latitudine geografica senza dover ricorrere ad alcun schemanumerico.

In realtà questa formula non è esatta ma approssima molto bene l’angoloϕ sotto certecondizioni che mostreremo tra poco.

Supponiamo ora di prendere un puntoappartenente alla superficie ellissoidica di coordinate

r,z individuate in un preciso piano meridiano che soddisfano quindi l’equazione (1.5.1) eandiamo a definire e individuare l’angoloθ per cui si ha

θ = arctanzarb

Ricordando il teorema (1.4.4) e applicando la trasformazioneb = a 1 − e2 ricaviamo larelazione

tanθ = 1 − e212 tanϕ (1.5.22)

Fin qui abbiamo compiuto semplici passaggi algebrici che all’apparenza non giustificanol’introduzione di questo nuovo angoloθ, anzi si potrebbe obbiettare persino che è perfettamenteinutile averlo introdotto, visto che in prossimità della superficie terrestre (cioè a quoteellissoidiche relativamente ‘basse’) per ricavare la latitudine geograficaϕ si potrebbe invertirepiù vantaggiosamente la relazione del teorema (1.4.4) anzichè la (1.5.22).

In realtà l’obbiettivo perseguito da Bowring era proprio quello di arrivare ad estrapolare unqualche tipo di formula per la latitudine geografica che fornisse dei valori sufficientementecorrettianche in regioni di spazio molto distanti dalla superficie terrestre.

Più avanti illustreremo graficamente l’andamento qualitativo degli errori commessi dallaformula di Bowring (1.4.5) e li confronteremo con quelli commessi applicando ingenuamentel’inversione della formula (1.4.4) decretando una volta per tutte la superiorità della prima sullaseconda.

Fatte queste basilari premesse possiamo finalmente indicare la strada ‘tortuosa’ (una dellepossibili) che portano alla scoperta della formula.

Dalla relazione (1.5.22) possiamo ricavare facilmente anche le funzioni goniometriche delseno e del coseno diθ:

sinθ = tanθ1 + tan2θ

=1 − e2 tanϕ

1 + 1 − e2 tan2ϕ

cosθ = 11 + tan2θ

= 11 + 1 − e2 tan2ϕ

Teniamo in serbo queste formule perchè ci serviranno tra pochissimo e andiamo a considerareadesso questa strana identità senza peraltro motivarne la sua comparsa :

tanϕ

1 + 1 − e2 tan2ϕ=

tanϕ + 1 − e2 tan3ϕ

1 + 1 − e2 tan2ϕ32

Come si vede , non c’è dubbio alcuno che si tratti effettivamente di una identità trigonometricain quanto il secondo membro è stato ottenuto moltiplicando edividendo banalmente il primomembro per la quantità 1+ 1 − e2 tan2ϕ.

In particolare il secondo membro si può scomporre efficaciemente nella somma :

tanϕ

1 + 1 − e2 tan2ϕ32

+1 − e2

32 tan3ϕ

1 − e2 1 + 1 − e2 tan2ϕ32

= cos3θ tanϕ + sin3θ1 − e2

dove si vede benissimo finalmente come abbiamo utilizzato proficuamente le formule del seno ecoseno diθ.

Quindi, l’identità introdotta in precedenza ad hoc si può riscrivere più compattamente come

cos3θ tanϕ + sin3θ1 − e2

=tanϕ

1 + 1 − e2 tan2ϕ

Ricordiamo ancora una volta , semmai ce ne fosse bisogno , chel’identità scritta sopra è validasolo se gli angoliθ eϕ identificano un medesimo punto appartenente alla superficie terrestre.

A questo punto entra la vera novità del ragionamento di Bowring : estendiamo la definizionedi θ a qualsiasi punto di coordinate generiches,z anche esterne alla superficie ellissoidica

θ = arctanzasb

(1.5.23)

e dopo aver moltiplicato entrambi i membri della identità per il fattor comune−ae2, possiamoapplicare diligentemente la formula (1.5.12) e sostituirepertanto il secondo membro con laquantitàz− stanϕ

− ae2 cos3θ tanϕ − ae2 sin3θ1 − e2

≃ z− stanϕ (1.5.24)

La relazione (1.5.24) non è più una identità proprio perchè abbiamo adottato la definizione(1.5.23); la confusione nasce dal fatto di aver sostituito il parametror (che identifica l’ascissa diun punto che sta sulla superficie e ha latitudineϕ) con la generica coordinatas in (1.5.23), sicchèla (1.5.24) non poteva continuare più ad essere una identità.

Non solo: questa volta l’altezzaznon è subordinata alla ascissasproprio perchè abbiamorotto il vincolo alla superficie.

La ragione della estensione è ovvia: come avremmo potuto infatti calcolare il parametrocorrettor (e la sua corrispondente ordinata) conoscendo soltanto le coordinates,ze i parametriellissoidicia,e senza ricorrere all’uso di qualche metodo iterativo ?

Assumiamo comunque che la (1.5.24) continui a essere una identità ben consapevoli del fattoche risolvendola inϕ si avrà un errore di una certa entità sul reale valore della vera latitudine.

L’azzardo è grande, ma se riusciremo a dimostrare che la formula di Bowring è affettanonostante ciò da errori molto modesti, allora anche la (1.5.24) avrà piena cittadinanza.

Allora continuiamo a procedere nella deduzione di tale formula e riscriviamo l’identità(1.5.24) spostando un po’ di termini e mettendo bene in evidenza a secondo membro la secondaeccentricitàe′2 = e21 − e2−1 e il semiasse minoreb = a 1 − e2

stanϕ − ae2 cos3θ tanϕ = z+ a 1 − e2 e2

1 − e2 sin3θ

poi, raccogliendo la tangente diϕ a primo membro, riscriviamo il tutto più compattamente

s− ae2 cos3θ tanϕ = z+ be′2 sin3θ

da cui isoliamo finalmente la tangente e ricaviamo così la formula di Bowring

tanϕ = z+ be′2 sin3θs− ae2 cos3θ

(1.5.25)

Senza effettuare uno studio particolareggiato sulla funzione errore intendiamo comunque dareuna importante raccomandazione che consiste nel non utilizzare la formula di Bowring pereccentricità molto prossime all’unità in quanto gli scostamenti potrebbero essere molto forti.

Nel modello che stiamo trattando, cioè quello terrestre, ilproblema ovviamente non sussistepoichè l’eccentricità è molto piccola.

Il grafico riportato nella prossima pagina è il frutto di unasimulazione numerica effettuata alcomputer che riporta le curve degli errori commessi utilizzando la formula di Bowring a diverseinclinazioni della latitudineϕ.

Più precisamente sulle ascisse è stata riportata la quota ellissoidicah espressa in kilometri

(fino ad un massimo di 100.000 Km.) e sulle ordinate vengono indicate ovviamente le deviazioniangolari in gradi tra la latitudine vera e quella calcolata secondo Bowring.

Fig. (1.5.26)

Dal grafico si deduce che l’errore massimo si ha in corrispondenza della latitudineϕ ≃ 45° acirca 15.000 km di quota ellissoidica e vale approssimativamente 5⋅ 10−7 gradi, cioè meno di 2millesimi di secondo d’arco e ciò giustifica ampiamente il successo della formula di Bowring.

Per angoli superiori ai 45° l’errore tende ad assottigliarsi nuovamente ed è la ragione per cuile curve ad inclinazioni superiori non sono state rappresentate.

Che cosa possiamo dire invece degli errori commessi ad altezze ellissoidiche negative ?Evidentemente va detto che la risposta ha più importanza in ambito matematico che pratico

in quanto le perforazioni della crosta terrestre possono raggiungere al massimo uno o forse duekilometri e comunque il punto più basso degli oceani si colloca ad una profondità uguale a circa11 km che rapportata al raggio del pianeta vale percentualmente solo lo 0,17%.

Si è scoperto comunque con un altra simulazione al computer che fino ai 6300 Km.(ladimensione tipica del raggio del pianeta) la formula di Bowring si comporta ancora molto benegenerando errori dell’ordine dei 10−3 gradi , mentre esiste una banda compresa tra i 6300 - 6700Km in cui gli errori oscillano fortemente di svariate decinedi gradi rendendo la formulacompletamente inutilizzabile per poi, superata la soglia superiore, stabilizzarsi ad un valorecostante (ad es. a 9° tale valore era poco più di 18°.

La spiegazione di un tale comportamento consiste nel fatto che evidentemente a quoteellissoidiche negative superiori ai 6300 Km. le latitudinicorrette andrebbero calcolate seguendoun criterio di logicità, ad esempio supponiamo che il punto rappresentativo di una quotaellissoidica negativa (oltre i - 12.700 Km.) sia esterno al pianeta, allora la latitudine correttaandrebbe cercata facendo si’ che il medesimo punto appaia aduna quota ellissoidica positiva daqualche altra località.

A questo punto è interessante scoprire mediante simulazione al computer quale andamentohanno le curve di errore attribuite alla semplice inversione della formula (1.4.4) :

ϕ = arctan zs− se2 (1.5.27)

Nella figura sottostante si può notare che questa volta l’ordine di grandezza degli errori èintorno al decimo di grado a distanze molto elevate (per esempio a quota 36.000 Km. dove sitrovano le orbite geostazionarie) ed è decisamente troppo grande, specialmente in riferimento adapplicazioni satellitari per scopi ad uso civile o militarein cui la precisione è naturalmente unfattore essenziale.

Fig. (1.5.28)

Inoltre le informazioni contenute nei due grafici rivelanoaltre differenze sostanziali, infatti lecurve di errore secondo Bowring hanno un picco massimo ad unacerta quota dopodichè tendonoa diminuire incrementando la quota, mentre le curve di errore di questo tipo cresconoasintoticamente per quote sempre maggiori fino a qualche valore limite.

Questo comportamento non deve stupirci più di tanto perchè se torniamo a considerare lafigura (1.5.1) e immaginiamo di allontare il puntoPh indefinitamente lungo la normaleellissoidica, allora l’angoloψ (centrato in O e congiungentePh) tenderà ad eguagliare lalatitudine geograficaϕ.

Si dimostra abbastanza agevolmente che l’errore asintotico eϕ ad una certa latitudineϕ èdato da

eϕ = ϕ − arctantanϕ1 − e2 (1.5.29)

Alberto Cella24/08/2009