2016-02-12 2 views
2

Je suis en train d'implémenter les calculs du livre Practical Astronomy with your Calculator or Spreadsheet. Jusqu'à présent mes calculs donnent exactement le même résultat que les exemples de calculs du livre. Toutefois, en arrivant au §39 "Calcul des corrections pour la parallaxe", je rencontre une différence que je ne comprends pas.Erreurs dans les calculs mathématiques pour calculer les corrections pour la parallaxe de l'objet céleste en Java

La tâche est décrite comme suit:

À titre d'exemple, calculons l'ascension droite apparente et la déclinaison de la Lune le 26 Février 1979, à 16 heures TU 45m quand on l'observe à partir d'un emplacement 60 les coordonnées géocentriques étaient α = 22h 35m 19s et δ = -7 ° 41 '13' ', et la parallaxe horizontale équatoriale de la Lune était 1 ° 01' 09 ' '.

Le livre décrit la séquence de calcul comme suit: enter image description here

Pourtant mon issue de l'étape 7 est -31,993415, mais livre dit -31,99 415. Si je fais la mathématiques de l'étape 7 avec les valeurs du livre sur une calculatrice le résultat est -31,99 415 aussi, donc mon résultat semble être juste et le livre est faux ....

Je pourrais vivre avec ça mais là est une différence à l'étape 10 aussi. Mon résultat est -8,570634, le résultat des livres est -8,538165, une différence plutôt grande. J'ai lu l'étape 10 encore et encore pour voir s'il y a une erreur dans mon code, mais je ne le vois pas. Comme jusqu'à présent mes calculs et les calculs de livres sont exactement les mêmes que je suis coincé. Est-ce que je fais quelque chose de mal (préféré), ou ne le livre faire une erreur (espérons qu'il n'y a pas plus ...)

code Java Ma pour cette fonction est la suivante:

static EquatorialCoordinate parallax(EquatorialCoordinate body, ObserverLocation observer, ZonedDateTime zdt, double P) { 
    double Hd = 15d * raha(body.α, zdt, observer.λ); 
    step("α", body.α); 
    step("δ", body.δ); 
    step("φ", observer.φ); 
    step("λ", observer.λ); 
    step("h", observer.h); 
    step("H", Hd); 

    double H = toRadians(Hd); 
    Parallax ρ = parallax(observer.φ, observer.h); 

    step("P", P); 
    P = toRadians(P); 

    double δ = toRadians(body.δ); 
    double r = 1d/sin(P); 
    step("r", r); 
    double ρsinφ = ρ.sin; 
    double ρcosφ = ρ.cos; 
    step("ρcosφ'",ρcosφ); 
    step("ρsinφ'",ρsinφ); 
    double Δ = atan((ρcosφ * sin(H))/((r * cos(δ)) - (ρcosφ * cos(H)))); 
    step("Δ", toDegrees(Δ)); 
    H += Δ; 
    step("H'", toDegrees(H)); 
    Δ = toDegrees(Δ); 
    double α$ = body.α - (Δ/15d); 
    step("α'", α$); 

    double divident = (r * sin(δ)) - ρsinφ; 
    double divisor = (r * cos(δ) * cos(H)) - ρcosφ; 
    double δ$ = atan(cos(H) * (divident/divisor)); 
    δ$ = toDegrees(δ$); 
    step("δ'", δ$); 
    return new EquatorialCoordinate(α$, δ$); 
} 

La fonction "step" fait un simple printf formaté. La sortie de ce programme est:

  • α 22,588611
  • δ -7,686944
  • φ 50,000000
  • λ -100,000000
  • h 60,000000
  • H -31,642500
  • P 1,019167
  • r 56.221228
  • ρcosφ '0.644060
  • ρsinφ' 0,762422
  • Δ -0,350915
  • H '-31,993414
  • α' 22,612005
  • δ '-8,570634

Le résultat δ' est -8 ° 34' 14.28 "au lieu de -8 ° 32 '17"

J'ai remplacé ma valeur calculée de H' par la valeur comptable pour voir si le livre contient une erreur reportée, mais même si je le fais, la valeur est erronée.

Ainsi ... ma grande question est, est-ce que ma mise en œuvre est fausse (et où, je ne peux pas la voir), ou si les calculs des livres étaient erronés. (Editer :) La classe est annotée avec strictfp, en utilisant java.util.StrictMath.

+3

vos calculs sont Supposant correct, cela est très probablement le résultat d'erreurs répétées dans 'double' calculs. Essayez d'utiliser un 'BigDecimal'. – Tunaki

+0

J'ai une question. Comment avez-vous réussi à taper ces symboles dans le code java? – SomeDude

+2

Je seconde cela: vous avez quelques calculs lourds ici et les erreurs de précision dues à l'utilisation de nombres à virgule flottante vont s'additionner. Utiliser 'BigDecimal' devrait aider. – Thomas

Répondre

0

Vous écrivez

H += Δ; 

Ce qui change la valeur de H.

Ensuite, vous écrivez

double δ$ = atan(cos(H) * (divident/divisor)); 

qui utilise la nouvelle version de H quand il doit utiliser la valeur ancienne.

0

Grâce à @svasa, j'ai trouvé que le diviseur de l'étape 10 devrait contenir H, pas H '. Le code est correct:

static EquatorialCoordinate parallax(EquatorialCoordinate body, ObserverLocation observer, ZonedDateTime zdt, double P) { 
    double H = toRadians(15d * raha(body.α, zdt, observer.λ)); 
    P = toRadians(P); 
    Parallax ρ = parallax(observer.φ, observer.h); 
    double δ = toRadians(body.δ); 
    double r = 1d/sin(P); 
    double Δ = atan((ρ.cosφ * sin(H))/((r * cos(δ)) - (ρ.cosφ * cos(H)))); 
    double H$ = H + Δ; 
    double α$ = body.α - (toDegrees(Δ)/15d); 
    double δ$ = toDegrees(atan(cos(H$) * ((r * sin(δ) - ρ.sinφ)/(r * cos(δ) * cos(H) - ρ.cosφ)))); 

    return new EquatorialCoordinate(α$, δ$); 
}