;
; DER VON CANTERS OPTIMIERTE WAGNER VII (HAMMER-WAGNER)
; =====================================================
;
; Kurzname:       Canters Wagner VII-Optimierung
; Originalname:   Optimised Hammer-Wagner
; Autor:          Karlheinz Wagner 1949, Frank Canters 2002
; Quelle:         Karlheinz Wagner, Kartographische Netzentwürfe, Leipzig 1949, S. 205ff.
;                 Frank Canters, small-scale Map-Projection Design, London 2002, S. 125, 184ff. 
; Richtung:       Direkt-Transformation
;
; Das Programm übernimmt die Koordinaten eines Punktes (x/y) und transformiert
; diese in einen Punkt (x'/y').
;
; Literatur:
; Wagner: Kartographische Netzentwürfe, Leipzig: Bibliographisches Institut 1949
; Fiala: Mathematische Kartographie, Berlin: Verlag Technik 1957
; Canters: Small-scale Map Projection Design, London: Taylor & Francis 2002

; (C) Rolf Böhm 2006

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

;
; Laufende Koordinaten
;
	_name	Canters~Wagner~VII~(Hammer-Wagner)-Optimierung
	_var	phi		; Geographische Breite
	_var	lambda		; Geographische Länge
	_var	alpha		; Polarkoordinaten Winkel
	_var	delta		; Polarkoordianten Distanz
	_var	psi		; Parameter (eine Art phi)
	_var	delta0		; Mittelpunktslänge
	_var	t0		; temporär
	_var	t1		; temporär
	_var	t2		; temporär
	_var	t3		; temporär
	_var	t4		; temporär
	_var	scale		; Maßstabszahl
	_var	lambda0		; Mittelmeridianlänge
	_var	phi-pol		; Parallelkreis der die Pollinienlänge bestimmt
	_var	phi-lim		; Grenzbreite bis zu der die Flächenverzerrung kontrolliert wird
	_var	zfv		; zulässige Flächenverzerrung in Prozent
	_var	lambda-		; Letzte umbezifferte Länge
	_var	m1		; Das m1 von Wagner
	_var	m2		; Das m2 von Wagner
	_var	n		; Das n von Wagner
	_var	k		; Das k von Wagner
	_var	p		; Das p von Wagner
	_var	C5.62		; Das 5.6229 im Wagner bei Entwurf 20c
	_var	C2.61		; Das 2.6162 im Wagner bei Entwurf 20c
	_var	------		; End of Symbol table

; x, y, x', y', Cx', Cy', °(, (°, pi, pi/2 etc. sind vordefinierte globale Konstanten
;
; Initialisierung
; ===============
;
	tstne	initial	077$
; Dialog
	pause	Hinweis:~Dieses~Programm~rechnet~eine~Vorwärtstransformation.\\Es~muss~mit~einer~direkt~arbeitenden~Projection~engine~abgearbeitet~werden.
	input	scale	Maßstabszahl
	clip	scale	1	1E12

	input	lambda0	Mittelpunktslänge~in~Grad
	clip	lambda0	-180	180

	mov	phi-pol	53.4560161
	mov	lambda-	125.586
	mov	phi-lim	60
	mov	zfv	0
	mov	p	1
	div	p	0.4709

; in Bogenmass
	mul	phi-pol	°(
	mul	phi-lim	°(
	mul	lambda-	°(
; Konstanten (Konfigurationsvariablen) berechnen
	mov	t1	zfv
	div	t1	100		; Prozent
	add	t1	1
	div	t1	2
	acos	t1
	div	t1	phi-lim
	mov	m2	t1

	mov	t1	zfv	
	div	t1	100
	add	t1	1
	mov	t2	phi-lim
	cos	t2
	mul	t1	t2
	acos	t1
	div	t1	phi-lim
	mov	m2	t1

	mov	t1	phi-pol
	sin	t1
	mov	t2	pi/2
	mul	t2	m2
	sin	t2
	div	t1	t2
	mov	m1	t1

	mov	t1	phi-pol
	div	t1	2
	sin	t1
	mul	t1	p
	mov	t2	lambda-
	div	t2	2
	sin	t2
	div	t1	t2
	root	t1	2
	mov	k	t1

	mov	t1	lambda-
	div	t1	pi
	mov	n	t1

	mov	t1	k
	mul	t1	2
	mov	t2	m1
	mul	t2	m2
	mul	t2	n
	root	t2	2
	div	t1	t2
	mov	C5.62	t1

	mov	t1	2
	mul	t2	k		; Root(n m1 m2) steht noch in t2 drin!
	div	t1	t2
	mov	C2.61	t1
; Programm ist initialisiert
	mov	initial	1
077$:

;
; SIMD-Laufbereich
; ================

;
; Eigentlicher Entwurf, dieser direkt!
; ------------------------------------
;
	mov	lambda	x		; Geographische Länge
	mov	phi	y		; Geographische Breite

	sub	lambda	lambda0	

	cmod	lambda	-180	180

; x, y haben eine Doppelbedeutung:
; - einmal die verdefinierten RTA-Input-Koordinaten, also eigentlich Phi und Lambda,
; - dann aber auch die ebenen Kartenkoordinaten, die errechnet werden

;
; Umrechnung in Bogenmaß
; ----------------------
;
	mul	phi	°(
	mul	lambda	°(

; Netzentwurf rechnen
; -------------------

; von Lambda, Phi nach X, Y 

	mov	psi	phi
	mul	psi	m2
	sin	psi
	mul	psi	m1
	asin	psi			; psi fertig

	mov	t2	lambda
	mul	t2	n
	cos	t2
	mov	delta	psi
	cos	delta
	mul	delta	t2
	acos	delta			; delta fertig

	mov	alpha	psi
	sin	alpha
	mov	t3	delta
	sin	t3
	div	alpha	t3
	acos	alpha			; alpha noch nicht ganz fertig ... das Vorzeichen fehlt noch ...
	errjump	out

	mov	t0	lambda		; ... es ist dasselbe ...
	sgn	t0			; ... Vorzeichen wie lambda ...
	mul	alpha	t0		; ... jetzt ist alpha fertig.

	mov	t4	delta
	div	t4	2
	sin	t4			; Das ist jetzt das sin(delta/2)

	mov	x	alpha	
	sin	x
	mul	x	t4
	mul	x	C5.62		; x fertig

	mov	y	alpha	
	cos	y
	mul	y	t4
	mul	y	C2.61		; y fertig

;
; Maßstab, Kartenmittelpunkt etc. einrechnen
; ------------------------------------------
;
	mul	x	Rx'		; Erdradius
	div	x	scale		; Kartenmaßstab
	add	x	Cx'

	mul	y	Ry'
	div	y	scale
	add	y	Cy'
;
; Schlussarbeiten
; ---------------
;
111$:
	mov	x'	x
	mov	y'	y
	exit
	_end