Berechnung des D-Criteriums nach Jopek (1993)

Hallo liebe Spezialisten für Bahnberechnungen,
ich versuche mit Python das D-Criterium nach Jopek (1993) zu berechnen.

Die Ergebnisse der 3. Komponente "((2*s.sin(I21/2.0)) ** 2)" sehen allerdings komisch aus (>1). Ich vermute da ist in der Berechnung (vermutlich innerhalb von I21) etwas faul.
Oder muss das so? Hat jemand eine Idee?

Hier ein Beispiel:
Objekt 1                   Objekt 2                  D²                                    1. Komponente        2. Komponente                3. Komponente                4. Komponente
"(134340) Pluto"   "(136108) Haumea"   "1.79190658767661"   "0.004053996241"   "0.000289849536591"   "1.77152789670401"   "0.0160348451950088"



Hier ein paar Paper zu dem Thema:
https://academic.oup.com/mnras/article/327/2/623/1041222
http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1982BAICz..33..104K&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf
http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1963SCoA....7..261S&data_type=PDF_HIGH&whole_paper=YES&type=PRINTER&filetype=.pdf
https://academic.oup.com/mnras/article/455/4/4329/1273482

import scipy as s
import numpy as n 

        q1 = a1 * (1.0 + e1)
        q2 = a2 * (1.0 + e2)

        I21 = s.arccos( s.cos(i1) * s.cos(i2) + s.sin(i1) * s.sin(i2) * s.cos(O2 - O1))
        if n.absolute(O2 - O1) > 180:
            II21 = w2 - w1 - 2*s.arcsin( s.cos((i2 + i1)/2.0) * s.sin((O2 - O1)/2.0) * 1.0/s.cos(I21/2.0) )
        else:
            II21 = w2 - w1 + 2*s.arcsin( s.cos((i2 + i1)/2.0) * s.sin((O2 - O1)/2.0) * 1.0/s.cos(I21/2.0) )
           
        D2 = ((e2-e1) ** 2) + (((q2 - q1) / (q2 + q1)) ** 2) + ((2*s.sin(I21/2.0)) ** 2) + (((e2+e1)/2.0) ** 2) * ((2*s.sin(II21/2.0)) ** 2)