;
; SPHÄRISCHE ERDKUGELTRANSFORMATION IN FREIE LAGE
; ===============================================
;
; Name:           Sphärische Kugeltransformation
; Richtung:       Inverse Transformation
;
; Das Gradnetz der Erde wird in eine andere Lage auf der Erde transformiert.
; Nach der Transformation befindet sich die Erdachse in beliebiger Lage, d. h.
; Nord- und Südpol können beliebig verändert werden.
;
; Es handelt sich um eine sog. Metatransformation, d. h. eine Vor-Transformation.
; Es wird nicht in eine Kartenebene projiziert, vielmehr bleibt das Modell
; Erdkugel in »geographischer Projektion« bestehen. Zielbild muss eine Erdkugel sein, d. h. ein
; Bild mit Geokoordinaten (-180, -90) ... (180, 90).
;
; Man kann das Programm insbesondere dazu verwenden, um Programme, die
; Netzentwürfe in normaler Lage rechnen, schiefe (freie) Lagen zu berechnen.
;
; Parameter für Briesemeister (mit Transformator gerechnet): -45, 10, 180.
; Parameter für Nordic Projection (mit Transformator gerechnet): -45, 0, 180.
; Parameter für Canters C8-4 (Canters 2002, S. 209f., Centre 61, 24, 0): -29, 24, 180.
; Parameter für Canters C8-5 (Canters 2002, S. 211f., Centre 45, 20, 0): -45, 20, 180. 
; Parameter für Canters C8-6 (Canters 2002, S. 212f., Pol 30, -140, 30 verdreht): -30, 40, 150.
;
; Das Programm übernimmt die Koordinaten eines Punktes (x/y) und transformiert
; diese in einen Punkt (x'/y').
;
; Erdachse und eine Längenverdrehung werden abgefragt, einen Maßstab gibt es nicht,
; die Erde bleibt die Erde. Die Zielbildgeometrie wird dem Sekundäroperanden
; oder einem Fixbild entnommen.
;
; Literatur:
; Wagner: Kartographische Netzentwürfe, Leipzig: Bibliographisches Institut 1949
; Fiala: Mathematische Kartographie, Berlin: Verlag Technik 1957
; Evenden: A Comprehensive Library of Cartographic Projection Functions 2004
;

; (C) ROlf Böhm 2004

; Benutzte Variablen
; ==================
;
; Die Variablennamen entsprechen weitgehend denen von Karlheinz Wagner,

;
; Laufende Koordinaten
;
	_name	Sphärische~Erdkugeltransformation
	_var	phi		; Geographische Breite
	_var	lambda		; Geographische Länge
	_var	delta		; Poldistanz/geographisch, auch Schiefdistanz
	_var	sinl		; Sinus(Lambda)
	_var	cosl		; Cosinus(Lambda)
	_var	sinp		; Sinus(Phi)
	_var	cosp		; Cosinus(Phi)
	_var	phi'		; Kartographische Breite
	_var	lambda'		; Kartographische Länge in 4 Varianten
	_var	lambda"
	_var	lambda-
	_var	lambda=
	_var	d11
	_var	d12
	_var	d21
	_var	d22
	_var	dx		; Differenzen »gekreuzt«
	_var	d=		; Differenzen »even«
	_var	t1		; Temporäre Variablen
	_var	t2
	_var	t3
	_var	t4
;
; Konstanten der Transformation
;
	_var	phiP		; Breite des Pols
	_var	lambdaP		; Länge des Pols
	_var	lambda0		; Geogr. Länge des Bildmittelpunktes
	_var	phiP(		; dito im Bogenmaß
	_var	lambdaP(	; dito im Bogenmaß
	_var	lambda0(	; dito im Bogenmaß

	_var	sinphiP		; Sin(phiP)
	_var	cosphiP		; cos(phiP)
	_var	-pi
	_var	2pi

000$:
;
; x, y, x', y', °(, (°, pi, pi/2 etc. sind vordefinierte globale Konstanten
;
; Initialisierung
; ===============
;
	tstne	initial	77$
; Dialog

	input	phiP	Geographische~Breite~der~neuen~Erdachse~in~Grad
	input	lambdaP	Geographische~Länge~der~neuen~Erdachse~in~Grad
	input	lambda0	Rotation~um~den~neuen~Pol~in~Längengraden
; Eingegebene Werte auf Min/Max bringen und in Bogenmaß umrechnen
	clip	phiP	-90	90
	clip	lambdaP	-180	180
	clip	lambda0	-180	180
	mov	phiP(	phiP
	mul	phiP(	°(
	mov	lambdaP(	lambdaP
	mul	lambdaP(	°(
	mov	lambda0(	lambda0
	mul	lambda0(	°(
	mov	sinphiP	phiP(
	sin	sinphiP
	mov	cosphiP	phiP(
	cos	cosphiP
; Pi-Konstanten
	mov	2pi	pi
	add	2pi	pi
	mov	-pi	pi
	neg	-pi	pi
; Programm ist initialisiert
	mov	initial	1
77$:
;
; SIMD-Laufbereich
; ================

;
; Maßstab, Kartenmittelpunkt etc. einrechnen
; ------------------------------------------

; Weil es eine Metatransformation, die Koordinaten nicht so ...
;       sub     x       Cx'                ; Bildmittelpunkt
;       div     x       Rx'                ; Erdradius
;       mul     x       scale              ; Kartenmaßstab
;       sub     y       Cy'
;       div     y       Ry'
;       mul     y       scale
; ... sondern so vorbehandeln ...
	mul	x	°(		; Bogenmaß
	mul	y	°(		; Bogenmaß
;
; Eigentlicher Entwurf, dieser invers
; -----------------------------------
;
	mov	lambda	x		; Geographische Länge
	mov	phi	y		; Geographische Breite
;
; Zylinderverdrehung
; ------------------
;
	add	lambda	lambda0(		; im Bogenmaß!
;
; Umrechnung der schief (gedachten) phi/lambda in polständige phi'/lambda'-Koordinatenpaare
; -----------------------------------------------------------------------------------------
;
; Das Problem ist, dass die Rechnung über arcsin läuft und dieser doppeldeutig ist.
; Die Lösung: Die Sache wird ein zweitesmal mit einer Formal angegangen, die einen
; arccos liefert. Dieser ist nun zwar auch doppeldeutig. Es ist aber nun diejenige Lösung
; richtig, die sowohl arcsin als auch arccos liefern. Wegen Rundungsunsicherheiten wird
; nicht auf Identität, sondern auf Differenzminimum getestet.

; Vorab: Delta-Lambda, sin/cos von Delta-Lambda, Lambda und Phi berechnen

	mov	sinl	lambda
	sin	sinl		; Sin(Lambda)
	mov	cosl	lambda
	cos	cosl		; Cos(Lambda)
	mov	sinp	phi
	sin	sinp		; Sin(Phi)
	mov	cosp	phi
	cos	cosp		; Cos(Phi)

; phi' berechnen

	mov	t1	sinp		; phi' berechnen (Formeln Fiala, S. 79f.)
	mul	t1	sinphiP
	mov	t2	cosp
	mul	t2	cosphiP
	mul	t2	cosl
	add	t1	t2
	mov	phi'	t1
	asin	phi'

; Lambda (1): lambda' berechnen: Mittels Arcsin

	mov	t1	cosp
	neg	t1
	mul	t1	sinl
	mov	t2	phi'
	cos	t2
	div	t1	t2
	mov	lambda'	t1
	asin	lambda'

; Lambda (2): lambda- berechnen: Die zweite Lösung, die der Arcsin liefert

	tstlt	lambda'	201$
	mov	lambda-	pi
	sub	lambda-	lambda'
	jump	203$
201$:
	clr	lambda-
	sub	lambda-	pi
	sub	lambda-	lambda'
203$:

; Lambda (3): lambda" berechnen: Mittels Arccos

	mov	t1	cosphiP
	neg	t1
	mul	t1	sinp
	mov	t2	sinphiP
	mul	t2	cosp
	mul	t2	cosl
	mov	t3	phi'
	cos	t3
	add	t1	t2
	div	t1	t3
	mov	lambda"	t1
	acos	lambda"

; Lambda (4): lambda= berechnen: Die zweite Lösung, die der arccos liefert

	mov	lambda=	lambda"		; Gleichheitszeichen lies als Doppelstrich quer
	neg	lambda=

; Lambda (5): Ein Lambda aus demjenigen Paar ist richtig, welches die kleinste Differenz aufweist

; Differenzen für 2. Entscheidungsstufe berechnen
	mov	d11	lambda'
	sub	d11	lambda"
	abs	d11
	mov	d12	lambda'
	sub	d12	lambda=
	abs	d12
	mov	d21	lambda-
	sub	d21	lambda"
	abs	d21
	mov	d22	lambda-
	sub	d22	lambda=
	abs	d22
; Differenzen für 1. Entscheidungsstufe berechnen
	mov	dx	d12		; "gekreuzt"
	mul	dx	d21
	mov	d=	d11		; "even"
	mul	d=	d22
; 1. Entscheidungsstufe: gekreuzt oder even
	cmpgt	dx	d=	=richtig
; 2. Entscheidungsstufe: Einzellösung
xrichtig:
	cmpgt	d12	d21	21richtig
	jump	12richtig
=richtig:
	cmpgt	d11	d22	22richtig
	jump	11richtig
; Lösungen zusammensuchen
21richtig:
	mov	lambda'	lambda-
	jump	ok
22richtig:
	mov	lambda'	lambda-
12richtig:
11richtig:
ok:
;
; Polverdrehung
; -------------
;
	add	lambda'	lambdaP(		; im Bogenmaß!
	cmpgt	lambda'	-pi	10$
	add	lambda'	2pi
10$:
	cmpgt	lambda'	-pi	30$
	add	lambda'	2pi
30$:
	cmplt	lambda'	pi	40$
	sub	lambda'	2pi
40$:
	cmplt	lambda'	pi	50$
	sub	lambda'	2pi
50$:
;
; In Gradmaß umrechnen und Ausserhalbtest
; ---------------------------------------

	mul	phi'	(°
	mul	lambda'	(°
;       cmplt   lambda'        -180        out
;       cmpgt   lambda'        180        out
;
; Schlussarbeiten
; ---------------
;
	neg	phi'		; Aus irgendeinem Grund negativ
	mov	x'	lambda'
	mov	y'	phi'
	exit
	end
out:
	mov	x'	-9999
	mov	y'	-9999
	exit
	_end