;
; CANTERS OPTIMIERUNG DES WINKELTREUEN KEGELENTWURFES FÜR AMERIKA (P06)
; =====================================================================
;
; Nr. bei Wagner: 5
; Name:           Winkeltreuer Schnittkegelentwurf nach Lambert mit zwei längentreuen Parallelkreisen;
;                 in schiefachsiger Lage, winkeltreuer Schnittkegelentwurf (I),
;                 Canters Optimierung für beide Amerika (Canters Part. Projection P06)
; Variante:       Zentrierend
; Autor:          Lambert 1772, Canters 2002
; Richtung:       Inverse Transformation
;
; Winkeltreuer Schnittkegelentwurf I =
; Schnittkegelentwurf mit zwei längentreu abgebildeten Parallelkreisen =
; Hauptentwurf der winkeltreuen Schnittkegel. Wagner 1949 S. 65

; Polversetzend/(Selbst-)Zentrierend: Bei Kegelentwürfen ist das Rechenzentrum (der Pol) meist 
; nicht der Kartenmittelpunkt. Dies ist ein (selbst)zentrierendes Programm, welches den Pol
; so versetzt, dass der Mitte der Quellkarte zur Mitte der Zielkarte wird. Das Gegenstück bilden
; polversetzende Polgramme, bei denen der Polversatz "zu Fuss" abgefragt wird. 
;
; Anmerkung 1: Die Berührungsparallel oder die Schnittparallele haben  k e i n e n  
;              Einfluss auf die Zentrierung. Entscheidend ist der Quellbildmittelparallel
;
; Anmerkung 2: Auf die X-Lagerung hat die Selbstzentrierung ebenfalls keinen Einfluss, 
;              dies ist, wie bei allen anderen Programmen, mit lambda0 einzustellen.
;
; Anmerkung 3: Die Selbstzentrierung ist ganz schön aufwändig. Immerhin muss
;              dafür zusätzlich zu der Inversformel eine Vorwärtsformel implementiert
;              werden (wenngleich diese auch nur ein einziges Mal gerechnet wird)
;
; The Canters Parameters are:
;
;	Canters (Tab. 5.21):    Here implemented:
;	PhiP     29.35		     Phi0      29.35
;       LambdaP  -7.61               Lambda0  180.00 + (-7.61) = 171.24
;	Phi1     10.36		     Phi1      10.36
;	Phi2     34.16               Phi2      34.16
;       Meassured in Fig. 5.32       Lambda00  85.80 Empiric Value: Rotates Cape Dame Marie (The Haiti South West Cape) into the Map Centre
;
; Literatur:
; Wagner: Kartographische Netzentwürfe, Leipzig: Bibliographisches Institut 1949
; Fiala: Mathematische Kartographie, Berlin: Verlag Technik 1957
; Canters: Small-scale Map Projection Design. Lndon: Taylor&Fracis 2002.

; (C) Rolf Böhm Bad Schandau 2004, 2008


; Benutzte Variablen
; ==================

;  Die Variablennamen entsprechen weitgehend denen von Karlheinz Wagner,

	_name	Canters~Optimierung~des~winkeltreuen~Kegelentwurfes~für~Amerika~(P06)
	_var	phi		; Geographische Breite
	_var	lambda		; Geographische Länge
	_var	alpha		; Azimut ebene Polarkoord., auch Schiefazimut
	_var	delta		; Poldistanz/geographisch, auch Schiefdistanz
	_var	m		; Radius ebene Polarkoordinaten
	_var	n		; Winkelreduktionsfaktor
	_var	phi1		; 1. Schnittparallel
	_var	phi2		; 2. Schnittparallel
	_var	c		; Maßstabskonstante
	_var	phi0		; Mittlerer Parallel mit dem gerechnet wird
	_var	delta0		; Poldistanz des Berührungsparallels
	_var	delta1		; Poldistanz 1. Parallel
	_var	delta2		; Poldistanz 2. Parallel
	_var	cosd0		; Cosinus von Delta0
	_var	sigma		; Hemisphärenvorzeichen
	_var	T1		; Die dicke Konstante vor der Formel auf S. 63 oben
	_var	t1		; temporäre Variablen
	_var	t2
	_var	t3
	_var	t4
	_var	lambda0		; Geogr. Breite des Mittelmeridians
	_var	scale		; Kartenmaßstabszahl (also 1000000, nicht 1/1000000)
	_var	2pi		; 360 Grad
	_var	--		; Variable für Testausgaben
;
; Für Transformation in freie Lage
;
	_var	phi0
	_var	sinphi0
	_var	cosphi0
	_var	phi'				; Kartographische Breite
	_var	lambda'				; Kartographische Länge in 4 Varianten
	_var	lambda"
	_var	lambda-
	_var	lambda=
	_var	d11				; Differenzen
	_var	d12
	_var	d21
	_var	d22
	_var	dx				; Differenzen »gekreuzt«
	_var	d=				; Differenzen »even«
;
; x, y, x', y', Cx', Cy', Rx', Ry', °(, (°, pi, pi/2 etc. sind vordefinierte globale Variablen
;

; Initialisierung
; ===============
;
	tstne	initial	$077
; Dialog
	input	scale	Maßstab
	mov	phi1	10.36		; Canters
	mov	phi2	34.16		; Canters	
	mov	phi0	29.35		; Canters
	mov	lambda0	180
	sub	lambda0	7.61		; Canters Value. Rechenwert 180 + (-8.76) = 172.39
	mov	lambda00	85.80	; Das dreht Amerika in auf die Kartenmittelachse. This rotates Cape Dame Marie into Map Centre
; Eingegebene Werte auf Min/Max bringen
	clip	scale	1	1E12
	clip	lambda0	-180	180
	mul	lambda00	°(		; die wird im Bogenmaß gerechnet
; Sigma ermitteln
	mov	sigma	phi1
	add	sigma	phi2
	sgn	sigma
	mul	phi1	sigma
	mul	phi2	sigma
; SinCos von Phi0 bilden
	mul	phi0	°(
	mov	sinphi0	phi0
	sin	sinphi0
	mov	cosphi0	phi0
	cos	cosphi0
; Konstanten berechnen
	mul	phi1	°(
	mul	phi2	°(
	mov	delta1	pi/2
	sub	delta1	phi1		; delta1 fertig
	mov	delta2	pi/2
	sub	delta2	phi2		; delta2 fertig
; n = Cos(Delta0) und n berechnen (Formel (5), Wagner, S. 64)
	mov	t1	delta2
	sin	t1
	log	t1
	mov	t2	delta1
	sin	t2
	log	t2
	mov	t3	delta2
	div	t3	2
	tan	t3
	log	t3
	mov	t4	delta1
	div	t4	2
	tan	t4
	log	t4
	sub	t1	t2
	sub	t3	t4
	div	t1	t3
	mov	n	t1		; n = svw. Cos(Delta0)
	acos	t1			; Und so
	mov	delta0	t1		; wirds Delta0
; Konstante c berechnen (Formel 3. Wagner S. 63)
	mov	r1	delta1
	sin	r1

	mov	r3	delta1
	div	r3	2
	tan	r3
	power	r3	n
	mul	r3	n

	div	r1	r3
	mov	c	r1
; Pol-y berechnen (Für das Zentrieren)
	mov	r1	Cy		; "Phi" von Cy abholen. "Lambda" alias Cx sei immer 0 ...
	add	r1	29.35		; Don't forget that
	mul	r1	sigma		; ;;;;;
	mul	r1	°(		; "Phi" in Bogenmaß umrechnen
;
	mov	delta	pi/2		; Eigentlicher Entwurf: Phi (in r1) im m (alias Pol-y) umrechnen
	sub	delta	r1		; pi/2 - phi

	mov	m	delta	
	div	m	2
	tan	m
	power	m	n		; svw. cos(delta0)
	mul	m	c
	mov	Pol-y	m
;					; m in y umrechnen. Aus x=0 folgt, da m=sqrt(x²+y²) m=y!
	div	Pol-y	scale		; Maßstab einrechnen
	mul	Pol-y	Ry'             ; Erdradius einrechnen
	mul	Pol-y	sigma		; hemisphärenabhängig
	proof	Pol-y
; 2pi
	mov	2pi	pi
	mul	2pi	2
; Programm ist initialisiert
	mov	initial	1
$077:
;
; SIMD-Laufbereich
; ================

;
; Maßstab, Kartenmittelpunkt etc. einrechnen
; ------------------------------------------
;
	mov	r0	x		; Change x<-->y
	mov	x	y
	mov	y	r0

	sub	x	Cy'		; Bildmittelpunkt
	neg	x
	div	x	Rx'		; Erdradius
	mul	x	scale		; Kartenmaßstab

	sub	y	Cx'		; Bildmittelpunkt
	sub	y	Pol-y		; Polverschiebung
	div	y	Ry'		; Erdradius
	mul	y	scale		; Kartenmaßstab

	mul	x	sigma
	mul	y	sigma
;
; Kartesische Koordinaten x/y in Polarkoordinaten m/alpha umwandeln
; -----------------------------------------------------------------

	mov	alpha	x		; kartesische Koord. x y in Polarkoordinaten m alpha
	div	alpha	y
	atan	alpha
	tstgt	y	$151		; Der atan ist doppeldeutig!: Wenn y negativ dann ...
	add	alpha	pi		; 180° addieren
$151:
	neg	alpha		; Jetzt ist alpha berechnet, liegt aber noch falsch
	add	alpha	pi
	cmplt	alpha	pi	$153
	sub	alpha	2pi
$153:
	power	x	2
	power	y	2
	clr	m
	add	m	x
	add	m	y
	root	m	2
;
; Eigentlicher Entwurf, dieser invers
; -----------------------------------
;
	mov	lambda	alpha		; Geographische Länge
	div	lambda	n
					; Geographische Breite
	mov	delta	m
	div	delta	c		
	root	delta	n
	errjump	out
	atan	delta
	mul	delta	2
	mov	phi	pi/2
	sub	phi	delta
;
; Kontern, Post-Transformation Rotation
; -------------------------------------
;
	mul	lambda	sigma		; Südhemisphäre:
	mul	phi	sigma		; ggf. kontern
	add	lambda	lambda00	; Post-Transformation, in the inverse code a pre-transformation (!)	
	cmod	lambda	-180	180
;
; Transformation in schiefe (freie) Lage: Umrechnung der polständigen phi/lambda in schiefe 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: sin/cos 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	sinphi0
	mov	t2	cosp
	mul	t2	cosphi0
	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	cosphi0
	neg	t1
	mul	t1	sinp
	mov	t2	sinphi0
	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:
;
; In Gradmaß umrechnen und Ausserhalbtest
; ---------------------------------------
;
	neg	lambda'			; ** Empirisch bemerkt, ohne den neg wirds seitenverkehrt.
	mul	lambda'	(°
	mul	phi'	(°
;
; Schlussarbeiten
; ---------------
;
	mov	x'	lambda'
	mov	y'	phi'
; Lambda kreisen lassen
	add	x'	lambda0
	cmod	x'	-180	180
	exit
out:
	mov	x'	-9999
	mov	y'	-9999
	exit
	_end