Algoritmo di Battin per risolvere il problema di Lambert

Si procede per successive approssimazioni della grandezza "y", uguale al rapporto tra il settore O-P1-P2 e il triangolo sotteso O-P1-P2. Ecco nell'ordine ed in sintesi, i vari "step" del procedimento iterativo.

Sono noti: r1(x1,y1,z1), r2(x2,y2,z2), , dt = t2-t1
  1. Calcolo dei parametri ausiliari: ( c, s, , , L, r0p, m)

    c= (r1² + r2² - 2 · r1 · r2 · cos )½   (corda dell'arco = distanza P1-P2)
    s= (r1 + r2 + c) / 2   (semi-perimetro del triangolo di Lambert)
    = (r1 · r2)½ / s · cos ( / 2)   (parametro adimensionale)
    cos = (r1 · r2)½ · cos ( / 2) / [½ · (r1+r2)]   da cui,   L= tan ² ( / 2)
    r0p = ¼ · (r1 + r2 + 2 · (r1· r2)½ · cos ( / 2))
    m = · dt ² / (8 · r0p ³)   con = parametro gravitazionale = 1 in unità canoniche

    "L" e "m" sono i parametri di Gauss; "r0p" è il raggio nel punto medio dell'arco di parabola che unisce P1 a P2.

  2. Scelta del valore iniziale x0 dell'iterazione (per es. x0=0, oppure x0 = L), essendo il parametro x compreso tra [-1 e +infinito] {saltare al punto 4}

  3. Calcolo di x dopo la prima iterazione x= [(1 - L) / 2)² + m / y²]½ - (1 + L) / 2

  4. Calcolo della Funzione Ipergeometrica :   posto = x / [1 + (1 + x)½] ²

    gs3= 169/25/27· /(1 + 196/27/29· / (1 + 225/29/31 · / (1 + 0)))

    gs2= 100/19/21 · / (1 + 121/21/23 · / (1 + 144/23/25 · / (1 + gs3)))

    gs1= 49/13/15 · / (1 + 64/15/17· / (1 + 81/17/19 · / (1 + gs2)))

    gs0= 16/7/9· / (1 + 25/9/11· / (1 + 36/11/13· / (1 + gs1)))

    si ottiene

    = 8·(1 + sqrt(1 + x))/(3 + 1/(5 + + 9/7· /(1+gs0)))

  5. Calcolo dei parametri (h1,h2):   posto hden= (1 + 2 · x + L) · (4 · x + · (3 + x)), si ha:

    h1= (L + x)² · (1 + 3 · x + ) / hden

    h2= m · (x - L + ) / hden

  6. Calcolo della Funzione Ipergeometrica "ku":   posto B= 27 · h2 / (4 · (1 + h1)·(1 + h1)²)

    u= -B / 2 / (1 + (1 + B)½)

    fs3= 928/3591·u/(1 - 1054/4347·u/(1 - 266/1035·u/(1 - 0)))

    fs2= 418/1755·u/(1 - 598/2295·u/(1 - 700/2907·u/(1 - fs3)))

    fs1= 22/81·u/(1 - 208/891·u/(1 - 340/1287·u/(1 - fs2)))

    fs0= 4/27·u/(1 - 8/27·u/(1 - 2/9·u/(1 - fs1)))

    si ottiene

    ku= 1/3/(1 - fs0)

  7. Radice positiva "y" del polinomio di 3° grado di Gauss

    yiterativo= (1 + h1) / 3 · (2 + (1 + B)½ / (1 - 2 · u · ku²)

    Se | yi+1 - yi | <   (ad es., 1·10-6) fine del processo iterativo al punto 8), altrimenti saltare al punto 3)

  8. Elementi dell'orbita (a, p, e)

    Detto "yf" il valore di "y" dell'ultima iterazione, si ottiene:

    a= m · s · (1 + )² / (8 · x · yf²)

    p= 2 · r1 · r2 · (yf · (1 + x) · sin ( / 2 ))² / ((1 + )² · m · s)

    e= (1 - p / a)1/2

  9. Coefficienti lagrangiani (f, g, f', g')

    f= 1 - r2 · 2 · sin ²( / 2 ) / p

    g= r1 · r2 · sin ( ) / ( · p)½

    f'= (2 / p · sin ² ( / 2) - 1/r1 - 1/r2) · tan ( / 2) · ( / p)1/2

    g'= 1 - r1 · 2 · sin ² ( / 2) / p

  10. Componenti delle velocità ai due estremi dell'arco P1 e P2

    Vx1=(x2-f·x1)/g

    Vy1=(y2-f·y1)/g

    Vz1=(z2-f·z1)/g

    Vx2=f'·x1+g'·Vx1

    Vy2=f'·y1+g'·Vy1

    Vz2=f'·z1+g'·Vz1

    Moduli delle velocità:

    V1=(Vx1²+Vy1²+Vz1²)1/2

    V2=(Vx2²+Vy2²+Vz2²)1/2





© Cassiope@