;
; McBRYDE-THOMAS FLAT-POLAR QUARTIC
; =================================
;
; Name:           McBryde-Thomas Flat-Polar Quartic
; Proj4:	  +proj=mbtfpq
;                 Libproj4 (2005), S. 64, 5.2.27
; Quelle:         Gerald I. Evenden. Libproj4: A Comprehensive Library if Cartographic Projection Functions, March 2004
; Richtung:       Direkttransformation
;
; Maßstab und Berührungsparallel werden abgefragt.
; Die Zielbildgeometrie wird dem Sekundäroperanden oder einem Fixbild entnommen.
; Der Berührungspunkt wird bildmittig gesetzt.
;
; Auflösung algebraisch nicht möglich. Berechnung mit Newton-Raphson-Iteration
;
; Literatur:
; Wagner: Kartographische Netzentwürfe, Leipzig: Bibliographisches Institut 1949
; Fiala: Mathematische Kartographie, Berlin: Verlag Technik 1957

; (C) Rolf Böhm 2004

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

;
; Laufende Koordinaten
;
	_name	McBryde-Thomas~Flat-Polar~Quartic
	_var	phi		; Geographische Breite
	_var	lambda		; Geographische Länge
	_var	delta		; Poldistanz/geographisch, auch Schiefdistanz
	_var	psi		; Parameter Psi
	_var	oldphi		; The last phi	 
	_var	p		; Newton-Raphson-Iteration: 0th derivative
	_var	p'		; Newton-Raphson-Iteration: 1th derivative
	_var	v		; psi Increment (p/p') "Verbesserung"
	_var	i		; Iteration Counter 
;
; Konstanten der Transformation
;
	_var	scale		; Kartenmaßstabszahl (also 1000000, nicht 1/1000000)
	_var	lambda0		; Geogr. Länge des Bildmittelpunktes
	_var	%2		; root(2)
	_var	2%3		; 2*root(3)
	_var	1+%2/2		; 1+root(2)/2
	_var	2+%2^		; 2+root(2)^(-0.5)
	_var	3%2+6^		; 3*root(2)+6^(-0.5)
;
; 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
	input	lambda0	Mittelpunktslänge~in~Grad
; Eingegebene Werte auf Min/Max bringen
	clip	scale	1	1E12
	clip	lambda0	-180	180
; Konstanten
	mov	%2	2		; "%" = "root"
	root	%2	2

	mov	2%3	3
	root	2%3	2
	mul	2%3	2

	mov	1+%2/2	%2
	div	1+%2/2	2
	add	1+%2/2	1

	mov	2+%2^	%2		; "^" = ^(-0.5)
	add	2+%2^	2
	power	2+%2^	-0.5

	mov	3%2+6^	%2
	mul	3%2+6^	3
	add	3%2+6^	6
	power	3%2+6^	-0.5
	
; Programm ist initialisiert
	mov	initial	1
077$:

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

;
; Lambda und Phi übernehmen
; -------------------------
;
	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
; -------------------
;
; psi - Newton-Raphson-Iteration: While Increment v < 1e-12 and no timeout 
;
	cmpeq	oldphi	phi	1030$	; Only 1 x per line
	mov	psi	phi		; Start psi is phi
	clr	i			; Iteration counter
	mov	oldphi	phi	
1010$:
	inc	i			; New cycle 
	cmpgt	i	1000	1030$	; Timeout

	clr	p			; 0th derivative:

	mov	r0	psi
	div	r0	2
	sin	r0
	add	p	r0		; p = sin(psi/2) ...

	mov	r1	psi
	sin	r1
	add	p	r1		; + sin(psi) ...

	mov	r2	phi
	sin	r2
	mul	r2	1+%2/2
	sub	p	r2		; - (1+root(2)/2) * sin(phi)
;
	clr	p'			; 1th derivative

	mov	r3	psi
	div	r3	2
	cos	r3
	mul	r3	0.5
	add	p'	r3		; p' = 0.5 * cos(psi/2) ...

	mov	r4	psi
	cos	r4
	add	p'	r4		; + cos(psi)	
;		
	mov	v	p
	div	v	p'		; v = p/p'
;
	sub	psi	v		; a new -- better -- psi
;
	mov	r7	v		; Compute error
	abs	r7
	cmpgt	r7	1e-12	1010$	; If abs(v) lt 1e-12 new Iteration
1030$:
;
; x
;
	mov	r0	psi
	div	r0	2
	cos	r0

	mov	x	psi
	cos	x
	mul	x	2
	add	x	1
	div	x	r0
	mul	x	3%2+6^
	mul	x	lambda
;
; y
;
	mov	y	psi
	div	y	2
	sin	y
	mul	y	2%3
	mul	y	2+%2^
;
; 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
;	mov	x'	psi		; Test
;	mul	x'	(°		; Test	
;	mov	y'	i		; Test
	exit
	_end