;
; McBRYDE THOMAS FLAT POLAR SINUSOIDAL
; ====================================
;
; Name:           McBryde-Thomas Flat-Polar Siusoidal
; Proj4:	  +proj=mbtfts
;                 Libproj4, S. 54, 6.2.1, 3/3 Sinusopial
; Quelle:         Gerald I. Evenden. Libproj4: A Comprehensive Library if Cartographic Projection Functions, March 2004
; Richtung:       Inverse Transformation
;
; Das Programm übernimmt die Koordinaten eines Punktes (x/y) und transformiert
; diese in einen Punkt (x'/y').
;
; x/y sind ebene Zielpunktkoordinaten, x'/y' geben die geogr. Breite und Länge
; der Position auf der Quell-Erdkugel, auf der der Zielpunkt gelesen werden
; kann.
;
; Maßstab und Berührungspunkt werden abgefragt.
; Die Zielbildgeometrie wird dem Sekundäroperanden oder einem Fixbild entnommen.
; Der Berührungspunkt wird bildmittig gesetzt.
;
; Literatur:
; Wagner: Kartographische Netzentwürfe, Leipzig: Bibliographisches Institut 1949
; Gerald I. Evenden, USGS: Libproj4. A Comprehensive Library if Cartographic Projection Functions, March 2004
;
; Formelzeichen u. a. aus Wagner-Quellen

; (C) Rolf Böhm 2004

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

;
; Laufende Koordinaten
;
	_name	McBryde-Thomas~Flat-Polar~Sinusoidal
	_var	phi		; Geographische Breite
	_var	lambda		; Geographische Länge
	_var	theta		; the parameter in the Evenden
;
; Konstanten der Transformation
;
	_var	lambda0		; Geogr. Länge des Bildmittelpunktes
	_var	scale		; Kartenmaßstabszahl (also 1000000, nicht 1/1000000)
	_var	m		; see Evenden
	_var	m+1		; m+1
	_var	n		; see Evenden 	
	_var	C		; see Evenden. In Wagner its the n
	_var	t1		; temp
	_var	t2		; temp
;
; x, y, x', y', Cx', Cy', °(, (°, pi, pi/2 etc. sind vordefinierte globale Konstanten
;
; Initialisierung
; ===============
;
	tstne	initial	077$
; Dialog
	input	scale	McBryde-Thomas~Flat-Polar~Sinusoidal\\Maßstabszahl
	input	lambda0	Mittelpunktslänge~in~Grad
; Eingegebene Werte auf Min/Max bringen
	clip	scale	1	1E12
	clip	lambda0	-180	180
; Konstanten berechnen
	mov	m	0.5		; Eckert IV: 1
	mov	m+1	m
	add	m+1	1

	mov	n	pi/4		; Eckert IV: pi/2
	add	n	1

	mov	t1	pi		; Eckert IV: C = 2/root(2+pi)
	add	t1	4
	mov	C	6
	div	C	t1
	root	t1	2
; Programm ist initialisiert
	mov	initial	1
077$:
;
; SIMD-Laufbereich
; ================
;
; Maßstab, Kartenmittelpunkt etc. einrechnen
; ------------------------------------------

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

	sub	y	Cy'
	div	y	Ry'
	mul	y	scale
;
; Eigentlicher Entwurf, dieser invers
; -----------------------------------
;
; The same code as Eckert IV but other Constants n m C
;
	mov	theta	y
	div	theta	C		; parameter theta ok

	mov	t1	theta		; only in mcbryde-th: theta check. 
	abs	t1			; without it:
	cmpgt	t1	pi/2	out	; some »additional world images«		

	mov	t1	theta
	mul	t1	m
	mov	t2	theta
	sin	t2
	add	t1	t2
	div	t1	n
	asin	t1
	errjump	out
	mov	phi	t1
	
	mov	lambda	x		; Geographische Länge
	mul	lambda	m+1
	mov	t1	theta
	cos	t1
	add	t1	m
	mul	t1	C
	div	lambda	t1
;
; In Gradmaß umrechnen und Ausserhalbtest
; ---------------------------------------
	mul	phi	(°
	mul	lambda	(°
	mul	theta	(°
	cmplt	phi	-90	out
	cmpgt	phi	90	out
	cmplt	lambda	-180	out
	cmpgt	lambda	180	out

;
; Schlussarbeiten
; ---------------
;
	mov	x'	lambda
	add	x'	lambda0
	mov	y'	phi
	mov	z'	theta		; test


	cmpgt	x'	-180	10$
	add	x'	360
10$:
	cmpgt	x'	-180	30$
	add	x'	360
30$:
	cmplt	x'	180	40$
	sub	x'	360
40$:
	cmplt	x'	180	50$
	sub	x'	360
50$:
	
	exit
out:
	mov	x'	-9999
	mov	y'	-9999
	exit
	_end