; ; WORLD PROJECTION ON REGULAR POLYHEDRA: ICOSAHEDRON (POLAR) ; ========================================================== ; ; Nr. bei Wagner: 14b ; Name: World projection on regular Polyhedra: Icosahedron (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/Icosahedron~(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 th ; triangle heigth _var pir ; pentagon inradius _var ch ; common triangle heigth (0.5*root(3)) _var gr ; golden ratio, pentagons diagonale (0.5*(1+root(5)) _var mh ; "main heigth" in triangle _var zh ; "Z (central) heigth" in triangle _var beta1 ; icosahedrons lower faces inclination _var -beta1 _var beta2 ; icosahedrons upper faces inclination _var -beta2 _var zeta ; angle between 2 icosahedron faces _var tau1 ; small hexagon angle _var -tau _var tau2 ; large hexagon angle _var -tau2 _var hcr ; hexagon circumradius _var t1 ; hexagon azimuths I = lower faces longitudes - north _var t2 _var t3 _var t4 _var t5 _var t6 _var u1 ; hexagon azimuths II = lower faces longitudes - south _var u2 _var u3 _var u4 _var u5 _var u6 ; _var a ; Face triangle sind length _var h ; Face triangle heigth ; _var 5h/6 ; multiples and fractions _var h/6 _var -h/6 _var -5h/6 _var -9a/4 _var -7a/4 _var -5a/4 _var -3a/4 _var -a/4 _var a/4 _var 3a/4 _var 5a/4 _var 7a/4 _var 9a/4 _var -pi/2 _var pi/3 _var -pi/3 _var 2pi/3 _var -2pi/3 _var -4pi/5 _var -3pi/5 _var -2pi/5 _var -pi/5 _var pi/5 _var 2pi/5 _var 3pi/5 _var 4pi/5 _var 5pi/5 ; ; Cube face selection ; _var inr ; Inradius _var fir ; Pentagon inradius _var fcr ; Face circum radius _var snoff ; Solid number of faces _dim FCCL 20 ; Face center central lambda array _dim FCCP 20 ; Face center central phi array _dim FCCR 20 ; Rotation array (not used) _dim FCSX 20 ; Face center shift X (in inradius units) array _dim FCSY 20 ; 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 ; ; INITIALISATION CODE ; =================== ; tstne initial 077$ ; ; Constants I: First ; mov ch 3 root ch 2 div ch 2 ; ch = Common (60°-60°-60°-)triangle Heigth: (1/2)*root(3) ; mov gr 5 root gr 2 add gr 1 div gr 2 ; Golden Ratio: (1/2)*(1+root(5)) ; ; Constants II: Lower faces longitude ; ; General: There are 3 kinds of icosahedrons faces in polar (pole is a triangle centre) aspect: ; 2 "polar" faces (inclination is zero, trivial), ; 6 "upper" faces (polar faces neiughbours, simple) and ; 12 "lower" faces ... not simple ... to solve it, ; cut a icosahedron parallel to a polar face triangle in 6 lower face triangle ; centres B C E F H I. It gives a hexagon or triangle ADG with inscribed hexagon: ; ; D ; / \ ; C-N-E ; / \ ; / Z \ ; B F ; / \ / \ ; A - I--M--H - G ; ; AB, CD, DE, FG ... are 3 regular 60°-60°-60°-triangles, side length 1, also BI, CE, FH length. ; BC, EF, HI are pentagon diagonals, length golden ratio gr. ; AD, DG, GA length is (1+gr+1). ; ; Z is the triangle and hexagons central point. M is the AG side half. N is the CE side half. ; The builded triangles AMD, CND and IMZ are rectangular. ; mov r0 gr ; mh: triangle ADG "main heigth" add r0 2 mov mh ch mul mh r0 ; mh: main heigth = MD = ADG heigth ; mov zh mh ; zh: triangle AMZ heigth div zh 3 ; zh: MZ = "Z heigth" ; mov r1 mh ; tau1: 1st longitude in triangle NZE ... sub r1 ch : ... by Pothagoras ... sub r1 zh ; cathedus 1 = ZN = MD - ND - MZ. (ND is the "common heigth" ch) mov r7 r1 ; save power r1 2 mov r2 0.5 ; cathedus 2 = NE = 0.5 power r2 2 mov r0 r1 add r0 r2 root r0 2 ; hypothenusis length mov tau1 r7 ; cathedus 1 div tau1 r0 acos tau1 ; tau1 = NZE angle = "1st longitude angle" mov -tau1 tau1 neg -tau1 mov hcr r0 ; hypothenusis = hexagon circumradius (1st computation) ; mov r1 gr ; tau2: 2nd longitude in triangle MZI by Pythagoras ... div r1 2 ; cathedus 1 = IM = (1/2)*gr power r1 2 mov r2 zh ; cathedus 2 = MZ = "Z heigth" zh power r2 2 mov r0 r1 add r0 r2 root r0 2 ; hypothenusis length mov tau2 zh div tau2 r0 acos tau2 ; tau2 = MZI angle = "2nd longitude angle" mov -tau2 tau2 neg -tau2 mov hcr r0 ; hypopthenusis = hexagon circumradius (2nd computation) ; ; (Remark: tau1 + tau2 = 60°) ; ; Constants III: Lower faces latitude or inclination ; mov r0 5 root r0 2 add r0 1 mov r3 3 root r3 2 ; r3: root(3) mul r0 r3 div r0 4 ; r0: icosahedrons dual dodecahedron circumradius ; mov r1 gr div r1 2 ; r1: pentagons diagonale half ; div r1 r0 asin r1 mov beta1 pi/2 sub beta1 r1 sub beta1 r1 mov -beta1 beta1 neg -beta1 ; ; Constants IV: Upper faces longitude ; ; ... trivial ... -2pi/3, 0, 2pi/3 ... ; ; Constants V: Upper faces latitude or inclination ; mov zeta 5 ; zeta: polar faces inclination - upper faces inclination (formula see Wikipedia) root zeta 2 div zeta 3 neg zeta acos zeta ; zeta ; mov beta2 zeta sub beta2 pi/2 mov -beta2 beta2 neg -beta2 ; ; Constants VI: Tables tau1, tau2 multiples in two angle tables t1 ... t6 (north hemishpere), u1 ... u6 (south). ; clr r0 ; t1 ... t3 add r0 tau2 mov t1 r0 add r0 tau1 add r0 tau1 mov t2 r0 add r0 tau2 add r0 tau2 mov t3 r0 ; clr r0 ; u1 ... u3 add r0 tau1 mov u1 r0 add r0 tau2 add r0 tau2 mov u2 r0 add r0 tau1 add r0 tau1 mov u3 r0 ; mov t6 t1 ; t4 ... t6 neg t6 mov t5 t2 neg t5 mov t4 t3 neg t4 ; mov u6 u1 ; u4 ... u6 neg u6 mov u5 u2 neg u5 mov u4 u3 neg u4 ; ; Constants VI: Shift values ; ; Based in "horizontal" driection on triangle height hz and in "vertical" ; direction on triangle side a. Note: Now other unit definition, unit sphere now ; is icosahedrons insphere! ; mov r0 5 root r0 2 add r0 3 mov r3 3 root r3 2 mul r0 r3 mov a 12 div a r0 ; triangle side length a in icosahedrons inradius units ; mov r0 3 root r0 2 div r0 2 mov h a mul h r0 ; triangle heigth h in icosahedrons inradius units ; ; Constants VII: Fractions and multiples ; mov h/6 h ; h fractions and multiples div h/6 6 mov 5h/6 h/6 mul 5h/6 5 mov -h/6 h/6 neg -h/6 mov -5h/6 5h/6 neg -5h/6 ; mov r0 a ; a fractions and multiples div r0 4 mov -9a/4 r0 mul -9a/4 -9 mov -7a/4 r0 mul -7a/4 -7 mov -5a/4 r0 mul -5a/4 -5 mov -3a/4 r0 mul -3a/4 -3 mov -a/4 r0 mul -a/4 -1 mov a/4 r0 mul a/4 1 mov 3a/4 r0 mul 3a/4 3 mov 5a/4 r0 mul 5a/4 5 mov 7a/4 r0 mul 7a/4 7 mov 9a/4 r0 mul 9a/4 9 ; mov -pi/2 pi/2 ; pi fractions and multiples neg -pi/2 mov pi/3 pi div pi/3 3 mov 2pi/3 pi/3 mul 2pi/3 2 mov -pi/3 pi/3 neg -pi/3 mov -2pi/3 2pi/3 neg -2pi/3 mov pi/6 pi div pi/6 6 mov -pi/6 pi/6 neg -pi/6 ; ; Solid definition I: Face center lambda/phi on Earth table ; ; Polar faces have + - pi/2 latitudes and 0 or pi longitudes. ; Upper faces have + - beta2 latitudes and + - pi/3 multiples as longitude. ; Lower faces have + - beta1 latitudes and + - tau1 or tau2 multiples as longitudes. ; The tau1 tau2 multiples are tabled as t1 t2 t3 ... (North) and u1 u2 u3 ... (South). ; mov FCCL(1) t4 ; lower N mov FCCP(1) beta1 mov FCCL(2) -2pi/3 ; upper N mov FCCP(2) beta2 mov FCCL(3) 0 ; polar N North pole mov FCCP(3) pi/2 mov FCCL(4) 2pi/3 ; upper N mov FCCP(4) beta2 mov FCCL(5) t3 ; lower N mov FCCP(5) beta1 mov FCCL(6) u4 ; lower S mov FCCP(6) -beta1 mov FCCL(7) t5 ; lower N mov FCCP(7) beta1 mov FCCL(8) 0 ; upper N Greenwich meridian mov FCCP(8) beta2 mov FCCL(9) t2 ; lower N mov FCCP(9) beta1 mov FCCL(10) u3 ; lower S mov FCCP(10) -beta1 mov FCCL(11) u5 ; lower S mov FCCP(11) -beta1 mov FCCL(12) t6 ; lower N mov FCCP(12) beta1 mov FCCL(13) t1 ; lower N mov FCCP(13) beta1 mov FCCL(14) u2 ; lower S mov FCCP(14) -beta1 mov FCCL(15) pi ; upper S mov FCCP(15) -beta2 mov FCCL(16) -pi/3 ; upper S mov FCCP(16) -beta2 mov FCCL(17) u6 ; lower S mov FCCP(17) -beta1 mov FCCL(18) u1 ; lower S mov FCCP(18) -beta1 mov FCCL(19) pi/3 ; upper S mov FCCP(19) -beta2 mov FCCL(20) pi ; polar S South pole mov FCCP(20) -pi/2 ; mov FCCR(1) -tau1 mov FCCR(2) -pi/3 mov FCCR(3) 0 mov FCCR(4) pi/3 mov FCCR(5) tau1 mov FCCR(6) -tau1 mov FCCR(7) -tau2 mov FCCR(8) 0 mov FCCR(9) tau2 mov FCCR(10) tau1 mov FCCR(11) -tau2 mov FCCR(12) -tau1 mov FCCR(13) tau1 mov FCCR(14) tau2 mov FCCR(15) 0 mov FCCR(16) -pi/3 mov FCCR(17) -tau1 mov FCCR(18) tau1 mov FCCR(19) pi/3 mov FCCR(20) 0 ; ; Solid definition II: Face center X/Y shift in the map table (in inradius units) ; mov FCSX(1) -9a/4 mov FCSY(1) 5h/6 mov FCSX(2) -5a/4 mov FCSY(2) 5h/6 mov FCSX(3) -a/4 mov FCSY(3) 5h/6 mov FCSX(4) 3a/4 mov FCSY(4) 5h/6 mov FCSX(5) 7a/4 mov FCSY(5) 5h/6 mov FCSX(6) -9a/4 mov FCSY(6) h/6 mov FCSX(7) -5a/4 mov FCSY(7) h/6 mov FCSX(8) -a/4 mov FCSY(8) h/6 mov FCSX(9) 3a/4 mov FCSY(9) h/6 mov FCSX(10) 7a/4 mov FCSY(10) h/6 mov FCSX(11) -7a/4 mov FCSY(11) -h/6 mov FCSX(12) -3a/4 mov FCSY(12) -h/6 mov FCSX(13) a/4 mov FCSY(13) -h/6 mov FCSX(14) 5a/4 mov FCSY(14) -h/6 mov FCSX(15) 9a/4 mov FCSY(15) -h/6 mov FCSX(16) -7a/4 mov FCSY(16) -5h/6 mov FCSX(17) -3a/4 mov FCSY(17) -5h/6 mov FCSX(18) a/4 mov FCSY(18) -5h/6 mov FCSX(19) 5a/4 mov FCSY(19) -5h/6 mov FCSX(20) 9a/4 mov FCSY(20) -5h/6 ; ; Solid definition III: Main values ; mov snoff 20 ; solid number of faces. Octahedron: 8 mov inr 1 ; insphere is the unit sphere mov fcr h ; face circumradius mul fcr 2 div fcr 3 ; ; 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 ; iot runs without rotations ; 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