;
; WINKELTREUER SCHNITTKEGELENTWURF IN NORMALER (POLARER) LAGE NACH LAMBERT MIT ZWEI LÄNGENTREUEN PARALLELKREISEN
; ==============================================================================================================
;
; Nr. bei Wagner: 5
; Name:           Winkeltreuer Schnittkegelentwurf nach Lambert mit zwei längentreuen Parallelkreisen
;                 Winkeltreuer Schnittkegelentwurf (I)
; Variante:       Zentrierend
; Autor:          Lambert 1772
; Quelle:         Karlheinz Wagner, Kartographische Netzentwürfe, Leipzig 1949. Es ist dies der 
;                 zweite Wagnersche Schnittkegel von Seite 65
; Richtung:       Direkttransformation
;
; Literatur:
; Wagner: Kartographische Netzentwürfe, Leipzig: Bibliographisches Institut 1949
; Fiala: Mathematische Kartographie, Berlin: Verlag Technik 1957

; (C) Rolf Böhm Bad Schandau 2004

; Winkeltreuer Schnittkegelentwurf I:

; Schnittkegelentwurf mit zwei längentreu abgebildeten Parallelkreisen.

; Hauptentwurf der winkeltreuen Schnittkegel. Wagner 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. 
;
; 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 mit lambda0 einstellen.

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

;
; Laufende Koordinaten
;
	_name	Winkeltreuer~Schnittkegelentwurf~(zentrierend)
	_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 (Cos(Delta0))
;
; Konstanten der Transformation
;
	_var	scale		; Kartenmaßstabszahl (also 1000000, nicht 1/1000000)
	_var	phi1		; Geogr. Breite des 1. Schnittparallels (besser Grenzparallels)
	_var	phi2		; Geogr. Breite des 2. Schnittparallels (besser Grenzparallels)
	_var	lambda0		; Geogr. Breite des Mittelmeridians
;
	_var	Pol-y		; Polversatz nach Nord
	_var	delta1		; Poldistanz des 1. Parallels
	_var	delta2		; Poldistanz des 2. Parallels
	_var	delta0		; Poldistanz des synthetischen Mittelparallels
	_var	cosd0		; dessen Cosinus
	_var	C	
;
; x, y, x', y', Cx', Cy', Rx', Ry', °(, (°, 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	phi1	1.~Schnittparallelbreite~in~Grad
	input	phi2	2.~Schnittparallelbreite~in~Grad
	input	lambda0	Berührungspunktlänge~in~Grad
; Eingegebene Werte auf Min/Max bringen
	clip	scale	1	1E12
	clip	phi1	-89.9999	89.9999
	clip	phi2	-89.9999	89.9999
	clip	lambda0	-180	180
; Konstanten berechnen
	mul	phi1	°(
	mov	delta1	pi/2		; delta1
	sub	delta1	phi1		
;
	mul	phi2	°(
	mov	delta2	pi/2		; delta2
	sub	delta2	phi2		
;					  Nun n und delta0 nach Formel Wagner S. 64, Gl. (5)
	mov	r1	delta1
	sin	r1
	log	r1
	mov	r2	delta2
	sin	r2
	log	r2
	sub	r2	r1
	mov	r3	delta1
	div	r3	2
	tan	r3
	log	r3
	mov	r4	delta2
	div	r4	2
	tan	r4
	log	r4
	sub	r4	r3
	mov	cosd0	r2
	div	cosd0	r4		; cosd0
	mov	delta0	cosd0
	acos	delta0			; delta0
; Nun gewöhnlich Winkeltreu (nun mit dem errechneten delta0)
	mov	n	cosd0		; n 
; 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 für Selbstzentrierung berechnen.
	mov	r0	Cy		; "Phi" von Cy abholen. "Lambda" alias Cx sei immer 0 ...
	mul	r0	°(		; "Phi" in Bogenmaß umrechnen
;					; Delta berechnen
	mov	delta	pi/2
	sub	delta	r0
;					; Eigentlicher Entwurf delta nach m (lambda ist immer 0)
	mov	m	delta		; Radius m
	div	m	2
	tan	m
	power	m	n
	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
	proof	Pol-y			; Kontrollausgabe
; 2pi
	mov	2pi	pi
	mul	2pi	2
; Programm ist initialisiert
	mov	initial	1
$077:
;
; SIMD-Laufbereich
; ================
;
	mov	lambda	x
	sub	lambda	lambda0
	cmod	lambda	-180	180
	mov	phi	y
;
; Umrechnung in Bogenmaß
; ----------------------
;
	mul	phi	°(
	mul	lambda	°(
;
; Eigentlicher Entwurf (in Polarkoordinaten)
; ------------------------------------------
;
	mov	alpha	lambda		; Azimut alpha
	mul	alpha	n

	mov	delta	pi/2		; delta
	sub	delta	phi

	mov	m	delta		; Radius m
	div	m	2
	tan	m
	power	m	n
	mul	m	C
;
; Polarkoordinaten m/alpha in kartesische Koordinaten
; ---------------------------------------------------
;
	mov	x	alpha		; kartesische Koord. x y in Polarkoordinaten m alpha
	sin	x
	mul	x	m
	mov	y	alpha
	cos	y
	neg	y			; ja, minus Cosinus!
	mul	y	m
;
; Maßstab, Kartenmittelpunkt etc. einrechnen
; ------------------------------------------
;
	mul	x	Rx'
	div	x	scale
	add	x	Cx'

	mul	y	Ry'
	div	y	scale
	add	y	Cy'
	add	y	Pol-y
;
; Schlussarbeiten
; ---------------
;
	mov	x'	x
	mov	y'	y
	exit
	_end