;
; FLÄCHENTREUE AZIMUTALPROJEKTION IN BELIEBIGER LAGE (SNYDER & VOXLAND)
; =====================================================================
;
; Name:           Flächentreue Azimutalprojektion in freier Lage
; Quelle:         Snyder & Voyland An Album of Map Projections, USGS Prof. Paper 1453 (1989)
; Richtung:       Direkt-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.

; (C) Dr.-Ing. Rolf Böhm 2005

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

;
; Laufende Koordinaten
;
	_name	Flächentreue~Azimutalprojektion~(Snyder) 
	_var	phi		; Geographische Breite
	_var	lambda		; Geographische Länge
	_var	z
	_var	K
	_var	sinphi
	_var	cosphi
	_var	sinphi0
	_var	cosphi0
	_var	sindl		; sin(Delta(Lambda))
	_var	cosdl		, cos(Delta(Lambda))
;
; Eingabeparameter
;
	_var	lambda0		; Geogr. Länge des Bildmittelpunktes
	_var	phi0		; bei Voxland phi1
	_var	scale		; Kartenmaßstabszahl (also 1000000, nicht 1/1000000)
	_var	--------
;
; 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	phi0	Berührungspunktbreite~in~Grad
	input	lambda0	Berührungspunktlänge~in~Grad
; Eingegebene Werte auf Min/Max bringen
	clip	scale	1	1E12
	clip	phi0	-90	90
	clip	lambda0	-180	180
; In Bogenmaß umrechnen
	mul	phi0	°(
	mul 	lambda0	°(
	mov	cosphi0	phi0				; Cos(phi0) braucht man nur einmal rechnen
	cos	cosphi0
	mov	sinphi0	phi0				; Sin(phi0) braucht man nur einmal rechnen
	sin	sinphi0
; Programm ist initialisiert
	mov	initial	1
077$:

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

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

;
; Umrechnung in Bogenmaß
;
	mul	phi	°(
	mul	lambda	°(
;
; Netzentwurf rechnen
; -------------------
;
; Winkelfunktionen
;
	mov 	sinphi	phi
	sin	sinphi
	mov	cosphi	phi
	cos	cosphi
	mov	sindl	lambda
	sub	sindl	lambda0
	sin	sindl
	mov	cosdl	lambda
	sub	cosdl	lambda0
	cos	cosdl

; Formel (112) S. 226

	mov	r0	sinphi0
	mul	r0	sinphi
	mov	r1	cosphi0
	mul	r1	cosphi
	mul	r1	cosdl
	clr	z	
	add	z	r0
	add	z	r1
	tsteq	z	-1	$out	; reject the point
	
;	acos	z			; z ist the cosine

; Formula (126) S. 228

	mov	r0	z	
	add	r0	1
	mov	K	2
	div	K	r0
	root	K	2

; Formula (114) S. 226

	mov	x	sindl
	mul	x	cosphi
	mul	x	K

; Formel (115) S. 226

	mov	r0	cosphi0	
	mul	r0	sinphi
	mov	r1	sinphi0
	mul	r1	cosphi
	mul	r1	cosdl
	mov	y	r0
	sub	y	r1
	mul	y	K
;
; 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
; ---------------
;
	mov	x'	x
	mov	y'	y
	exit
out:
$out:
	mov	x'	-9999
	mov	y'	-9999
	exit
	_end