;
; WORLD PROJECTION ON REGULAR POLYHEDRA: DODECAHEDRON (POLAR)
; ===========================================================
;
; Nr. bei Wagner: 14b
; Name:           World projection on regular Polyhedra: Dodecahedron (polar)
;                 Gnomonic Azimuthal Projection.
; Author:         Rolf Böhm
; Direction:      Inverse
;
; (C) Rolf Böhm Bad Schandau 2004, 2011, 2012
;
; DECLARATIONS
; ============
;
; Co-ordinates
;
	_name	World~projection~on~regular~polyhedra/Dodecahedron~(polar)
	_var	phi			; latitide
	_var	lambda			; longitude
	_var	alpha			; plane polar azimuth, oblique azimuth
	_var	delta			; plane polar distance, oblique distance
	_var	m			; plane polar radius
	_var	phi'			; cartographic latitude
	_var	lambda'			; cartographic longitude, solution 1 and common
	_var	lambda"                 ; cartographic longitude, solution 2
;
; Constants
;
	_var	phi0			; (oblique) central point latitude
	_var	sinphi0			; its sine
	_var	cosphi0			; its cosine
	_var	lambda0			; (oblique) central point longitude
	_var	lambda1			; rotation angle
	_var	scale			; map scale denominator (e. g. 1000000, not 1/1000000)
;
	_var	gr			; golden ratio
	_var	beta			; dodecahedrons face inclination
	_var	-beta			; 
;
	_var	a			; Pentagons side length
	_var	g			; Pentagons diagonale
	_var	r			; Pentagons circumradius
	_var	h			; Pentagons heigth
	_var	f			; Pentagons "first" 
	_var	c			; Pentagons "commma"
	_var	s			; Face pentagon "shift"
;
	_var	-c			; vertical shifts
	_var	hhc
	_var	-chh
;
	_var	s/2			; horizontal shifts
	_var	3s/2
	_var	5s/2
	_var	-s/2
	_var	-3s/2
	_var	-5s/2
;
	_var	-4pi/5
	_var	-2pi/5
	_var	-pi/5
	_var	pi/5
	_var	2pi/5
	_var	4pi/5
	_var	-pi
	_var	-pi/2
;
; Cube face selection
;
	_var	inr			; Inradius
	_var	fir			; Pentagon inradius
	_var	fcr			; Face circum radius
	_var	snoff			; Solid number of faces
	_dim	FCCL	12		; Face center central lambda array
	_dim	FCCP	12		; Face center central phi array
	_dim	FCCR	12		; Rotation array (not used)
	_dim	FCSX	12		; Face center shift X (in inradius units) array
	_dim	FCSY	12		; Face center shift y (in inradius units) array
	_var	d			; Distance on face
	_var	i			; Index of face
;
; Miscelleanous
;
	_var	sinl			; sin(lambda)
	_var	cosl			; cosine(lambda)
	_var	sinp			; sine(phi)
	_var	cosp			; cosine(phi)
	_var	dlambda			; delta lambda (=lambda-lambda0)
	_var	sindl			; sine delta lambda
	_var	cosdl			; cosine delta lambda
	_var	return			; return target address
	_var	initial			; program init flag
;
; x, y, x', y', Cx', Rx', °(, (°, pi, pi/2 etc. are pre-defined global variables
;
; INITIALIZATION CODE
; ===================
;
	tstne	initial	077$
; 
; Constants I: Dodecahedrons face inclination
;
;
	mov	beta	5
	root	beta	2
	div	beta	5
	neg	beta
	acos	beta
	sub	beta	pi/2
;
	mov	-beta	beta
	neg	-beta
;
; Constants II: Pentagon lines (if dodecahedrons insphere radius unit = 1)
;
; Defintions: Think a plane pentagon A B C D E with base bottom side A-B, top D, and C E as "first" (ridge) points, 
; Z is the pentagons center,
; M ist the base center point (mean A-B), 
; Q is the left side center point (mean E-A),
;  then
; a is the pentagons side length, e. g. A-B (pentagons base),
; g (golden ratio) is the pentagon diagonale length, e. g. C-E,
; r is the circumcircle radius, e. g. Z-D, 
; h ("heigth")is the incircle radius, e. g. M-Z,
; f ("first" or ridge) is the y value of C or E (meassured from base A-B) 
; c ("comma") is the y difference Z - Q, the vertical shift mean E-A point -- pentagon centre Z,
; s (horizontal "shift") is (a+g)/2, the horizontal shift between neighboured pentagons in the plane dodecahedrons graph.
;
	mov	gr	5
	root	gr	2
	add	gr	1
	div	gr	2	; golden ratio
;
	mov	a	2	; enumerator 
	mov	r0	5
	root	r0	2
	mul	r0	11
	add	r0	25
	div	r0	10
	root	r0	2	; denominator
	div	a	r0	; pentagons side if doedecahedrons insphere radius = 1
;
	mov	g	a
	mul	g	gr	; pentagons diagonale
;
	mov	r	5
	root	r	2
	mul	r	10
	add	r	50
	root	r	2
	div	r	10	
	mul	r	a	; pentagons circumradius (note: dodecahedrons inradius = 1)
;
	mov	h	a
	div	h	2
	mov	r0	pi
	mul	r0	3
	div	r0	10
	tan	r0		; sin(54°)
	mul	h	r0	; Pentagons inradius = a/2*sin(54°)
;
	mov	f	pi
	mul	f	2
	div	f	5
	sin	f
	mul	f	a	; pentagon first (=ridge heigth) = a*sin(72°)
;
	mov	c	h
	mov	r0	f
	div	r0	2
	sub	c	r0	; pentagons "comma": y difference between half first heifth anf heigth
;
	mov	s	a
	add	s	g
	div	s	2	; pentagons shift to (horizontal) next pentagon
;
; Constants III: shift values
;	
	mov	-c	c
	neg	-c
	mov	hhc	h
	add	hhc	h
	sub	hhc	c
	mov	-chh	hhc
	neg	-chh
;
	mov	s/2	s
	div	s/2	2
	mov	3s/2	s/2
	mul	3s/2	3
	mov	3s/2	s/2
	mul	3s/2	3
	mov	5s/2	s/2
	mul	5s/2	5
	mov	-s/2	s/2
	neg	-s/2
	mov	-3s/2	3s/2
	neg	-3s/2
	mov	-5s/2	5s/2
	neg	-5s/2
;
; Constants IV: pi fractions
;
	mov	-pi	pi
	neg	-pi
	mov	-pi/2	pi/2
	neg	-pi/2
;
	mov	r0	pi
	div	r0	5
;
	mov	-4pi/5	r0
	mul	-4pi/5	-4
	mov	-3pi/5	r0
	mul	-3pi/5	-3
	mov	-2pi/5	r0
	mul	-2pi/5	-2
	mov	-pi/5	r0
	mul	-pi/5	-1
	mov	pi/5	r0
	mul	pi/5	1
	mov	2pi/5	r0
	mul	2pi/5	2
	mov	3pi/5	r0
	mul	3pi/5	3
	mov	4pi/5	r0
	mul	4pi/5	4
;
; Solid definition I: Face center lambda/phi on Earth table
;
	mov	FCCL(1)		0
	mov	FCCP(1) 	pi/2
	mov	FCCL(2)		4pi/5
	mov	FCCP(2)		beta
	mov	FCCL(3) 	-4pi/5
	mov	FCCP(3)		beta
	mov	FCCL(4)		-2pi/5
	mov	FCCP(4)		beta
	mov	FCCL(5)		2pi/5
	mov	FCCP(5) 	beta
	mov	FCCL(6)		pi
	mov	FCCP(6)		-beta
	mov	FCCL(7) 	0
	mov	FCCP(7)		beta
	mov	FCCL(8)		3pi/5
	mov	FCCP(8)		-beta
	mov	FCCL(9) 	-3pi/5
	mov	FCCP(9)		-beta
	mov	FCCL(10)	-pi/5
	mov	FCCP(10)	-beta
	mov	FCCL(11)	pi/5
	mov	FCCP(11)	-beta
	mov	FCCL(12)	0
	mov	FCCP(12)	-pi/2
;
	mov	FCCR(1)		0
	mov	FCCR(2)		pi/5
	mov	FCCR(3)		-pi/5
	mov	FCCR(4)		-pi/5
	mov	FCCR(5)		pi/5
	mov	FCCR(6)		0
	mov	FCCR(7)		0
	mov	FCCR(8)		pi/5
	mov	FCCR(9)		-pi/5
	mov	FCCR(10)	-pi/5
	mov	FCCR(11)	pi/5
	mov	FCCR(12)	pi

;
; Solid definition II: Face center X/Y shift in the map table (in inradius units)
;
	mov	FCSX(1)		-3s/2	
	mov	FCSY(1)		hhc
	mov	FCSX(2)		s/2
	mov	FCSY(2) 	hhc
	mov	FCSX(3)		5s/2
	mov	FCSY(3)		hhc
	mov	FCSX(4)		-5s/2
	mov	FCSY(4)		c
	mov	FCSX(5)		-s/2		
	mov	FCSY(5)		c
	mov	FCSX(6)		3s/2
	mov	FCSY(6) 	c
	mov	FCSX(7)		-3s/2
	mov	FCSY(7)		-c
	mov	FCSX(8)		s/2
	mov	FCSY(8)		-c
	mov	FCSX(9)		5s/2
	mov	FCSY(9)		-c
	mov	FCSX(10)	-5s/2
	mov	FCSY(10)	-chh
	mov	FCSX(11)	-s/2
	mov	FCSY(11)	-chh
	mov	FCSX(12)	3s/2
	mov	FCSY(12)	-chh
;	
; Solid definition III: Main values
;
	mov	snoff	12		; solid number of faces. Dodecahedron: 12
	mov	inr	1		; insphere is the unit sphere
	mov	fcr	r		; face circumradius
;	
; Dialog
;
	input	scale	Map~scale~(denominator)
	clip	scale	1	1E12
;
; Program setted as initialized
;
	mov	initial	1
077$:

;
; SIMD-CODE
; =========
;
; Scale target geometry (0.1-mm-pixels) into centered unit sphere
; ---------------------------------------------------------------
;
	sub	x	Cx'			; Image center
	div	x	Rx'			; Earth radius
	mul	x	scale			; Map scale

	sub	y	Cy'
	div	y	Ry'
	mul	y	scale
;
; x, y now on unit sphere
;
; Detect the solids face index
; ----------------------------
;	
	mov	r0	snoff			; solid number of faces = index
	mov	d	1E90			; current distance ("infinite")
l_cycle:
	mov	return	l_return
	jump	.euclid				; returns the distance actual point - face center in r1
l_return:
	cmpgt	r1	d	l_continue
	mov	d	r1			; new distance
	mov	i	r0			; new index
l_continue:
	dec	r0
	tstgt	r0	l_cycle
	jump	l_end
l_end:
	cmplt	d	fcr	l_noout
	mov	x'	-9999
	mov	y'	-9999
	exit			; outside
l_noout:
;
; now contains ...
; i ... the current face index
; d ... the distance actual point - face center
;
; Get the detected face lambda0/phi0/rotation from table
; ------------------------------------------------------
;
	get	lambda0	FCCL	i		; lambda0
	get	phi0	FCCP	i		; phi0
	get	lambda1	FCCR	i		; rotation
;
	mov	cosphi0	phi0			; phi0 angle functions	
	cos	cosphi0
	mov	sinphi0	phi0		
	sin	sinphi0
;
; Page the single faces in a local system
; ---------------------------------------
;
	get	r0	FCSX	i
	sub	x	r0
	get	r2	FCSY	i
	sub	y	r2
;
; Projection I: Cartesian x/y into Polar co-ordinates m/alpha
; -----------------------------------------------------------
;
	mov	alpha	x
	div	alpha	y
	err	.	d_zero			; d_ error handler begin
	jump	d_okok
d_zero:
	mov	alpha	pi/2
	tstgt	x	d_end
	neg	alpha
	jump	d_end	
d_okok:
	atan	alpha
	tstgt	y	d_end			; If negative y then ...
	add	alpha	pi			; add 180°
d_end:						; d_ error handler end
151$:
	power	x	2
	power	y	2
	clr	m
	add	m	x
	add	m	y
	root	m	2
;
; Projection II: (Inverse) gnomonic projection
; --------------------------------------------

;  delta ... pole distance (90º-latitude)

	mov	lambda	alpha			; Longitude
	add	lambda	lambda1			; Rotation

	mov	delta	m			; Latitude
	atan	delta
	mov	phi	pi/2
	sub	phi	delta
;
; Projection III: Pole centered phi/lambda pair into oblique delta/alpha co-ordinates
; -----------------------------------------------------------------------------------
;
; alpha = (here) oblique longitude
; delta = (here) oblique pole distance (90º-latitude)
;
; Formulas: Fiala, Kartographische Netzentwürfe, Prague, p. 79f.

; compute constant angle functions

	mov	sinl	lambda
	sin	sinl			; sine(lambda)
	mov	cosl	lambda
	cos	cosl			; cosine(lambda)
	mov	sinp	phi
	sin	sinp			; sine(phi)
	mov	cosp	phi
	cos	cosp			; cosine(phi)

; phi'

	mov	r1	sinp		
	mul	r1	sinphi0
	mov	r2	cosp
	mul	r2	cosphi0
	mul	r2	cosl
	add	r1	r2
	mov	phi'	r1
	asin	phi'

; lambda 1st solution: arcsine

	mov	r1	cosp
	neg	r1
	mul	r1	sinl
	mov	r2	phi'
	cos	r2
	div	r1	r2
	mov	lambda'	r1

; lambda 2nd solution: arccosine 

	mov	r1	cosphi0
	neg	r1
	mul	r1	sinp
	mov	r2	sinphi0
	mul	r2	cosp
	mul	r2	cosl
	mov	r3	phi'
	cos	r3
	add	r1	r2
	div	r1	r3
	mov	lambda"	r1

; lambda by arcsine with the arccosine sign (using the FAI/full angle inversion instruction)

	asin	lambda'	lambda"
	add	lambda'	lambda0
;
; Oblique co-ordinates from arc into degree
; -----------------------------------------
;
	mul	phi'	(°
	mul	lambda'	(°
;
; Final Computations
; ------------------
;
	mov	x'	lambda'
	mov	y'	phi'

	add	x'	180
	cmpgt	x'	-180	3$
	add	x'	360
3$:
	cmplt	x'	180	4$
	sub	x'	360
4$:
	exit
;
; SUBROUTINES
; ===========
;
; .euclid - euclidean distance
; ----------------------------
;
; Returns the euclidean distance between 2 Points
; P1(x, y) and P2(FCSX(i), FCSY(i)) in r1. 
; x and y is a global co-ordinate pair.
; FCSX and FCSY are 2 global tables with x/y values, 
; r0 is the table index i.
;
.euclid:
	get	r1	FCSX	r0
	sub	r1	x
	power	r1	2
	get	r2	FCSY	r0
	sub	r2	y
	power	r2	2
	add	r1	r2
	root	r1	2
	jump	return
;
	_end