1 /*************************************** 2 * Copyright 2011, 2012 GlobWeb contributors. 3 * 4 * This file is part of GlobWeb. 5 * 6 * GlobWeb is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU Lesser General Public License as published by 8 * the Free Software Foundation, version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * GlobWeb is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with GlobWeb. If not, see <http://www.gnu.org/licenses/>. 18 ***************************************/ 19 20 /** 21 * Module which contains all the maths stuff 22 */ 23 define(['./HEALPixTables','./Long'], function(HealPixTables,Long) { 24 25 /**************************************************************************************************************/ 26 27 var HALF_PI = 3.14159265/2; 28 29 var lonLat2ang = function(lon, lat) 30 { 31 if ( lon < 0 ) 32 lon += 360; 33 34 var phi = lon * Math.PI / 180.; 35 36 var theta = ( -lat + 90. ) * Math.PI / 180.; 37 return [phi, theta]; 38 } 39 40 /**************************************************************************************************************/ 41 42 /** Returns the remainder of the division {@code v1/v2}. 43 The result is non-negative. 44 @param v1 dividend; can be positive or negative 45 @param v2 divisor; must be positive 46 @return Remainder of the division; positive and smaller than {@code v2} */ 47 var fmodulo = function(v1, v2) 48 { 49 if (v1>=0.) 50 return (v1<v2) ? v1 : v1%v2; 51 var tmp=v1%v2+v2; 52 return (tmp==v2) ? 0. : tmp; 53 } 54 55 /**************************************************************************************************************/ 56 57 var spread_bits = function(v) 58 { 59 return (HealPixTables.utab[ v &0xff]) | ((HealPixTables.utab[(v>>> 8)&0xff])<<16) 60 | ((HealPixTables.utab[(v>>>16)&0xff])<<32) | ((HealPixTables.utab[(v>>>24)&0xff])<<48); 61 } 62 63 /**************************************************************************************************************/ 64 65 var xyf2nest = function(ix, iy, face_num, order) 66 { 67 return ((face_num)<<(2*order)) + 68 spread_bits(ix) + (spread_bits(iy)<<1); 69 } 70 71 /**************************************************************************************************************/ 72 73 var loc2pix = function(order, phi, theta) 74 { 75 var nside = Math.pow(2, order); 76 var z = Math.cos(theta); 77 var phi = phi; 78 79 var loc = { 80 phi: phi, 81 theta: theta, 82 z: z 83 } 84 if (Math.abs(z)>(9./10.)) 85 { 86 loc.sth = Math.sin(theta); 87 loc.have_sth=true; 88 } 89 90 var inv_halfpi = 2./Math.PI; 91 var tt = fmodulo((phi*inv_halfpi),4.0);// in [0,4) 92 93 var za = Math.abs(z); 94 if (za<=2./3.) // Equatorial region 95 { 96 var temp1 = nside*(0.5+tt); 97 var temp2 = nside*(z*0.75); 98 99 var jp = Long.fromNumber(temp1 - temp2); 100 var jm = Long.fromNumber(temp1 + temp2); 101 var ifp = jp.shiftRightUnsigned(order); 102 var ifm = jm.shiftRightUnsigned(order); 103 var face_num; 104 if ( ifp.equals(ifm) ) 105 { 106 face_num = ifp.or(Long.fromInt(4)); 107 } 108 else 109 { 110 if ( ifp.lessThan(ifm) ) 111 { 112 face_num = ifp; 113 } 114 else 115 { 116 face_num = ifm.add(Long.fromInt(8)); 117 } 118 } 119 120 var nSideMinusOne = Long.fromNumber(nside-1); 121 var ix = jm.and( nSideMinusOne ); 122 var iy = nSideMinusOne.subtract( jp.and(nSideMinusOne) ); 123 124 return xyf2nest(ix.toInt(),iy.toInt(),face_num.toInt(), order); 125 126 } 127 else // polar region, za > 2/3 128 { 129 var ntt = parseInt( Math.min( 3, parseInt(tt) ) ); 130 var tp = tt-ntt; 131 var tmp = ( (za < (9./10.)) || (!loc.have_sth) ) ? 132 nside*Math.sqrt(3*(1-za)) : 133 nside*loc.sth/Math.sqrt((1.+za)/3.); 134 135 var jp = Long.fromNumber(tp*tmp); 136 var jm = Long.fromNumber((1.0-tp)*tmp); 137 var lNside = Long.fromNumber(nside); 138 var nSideMinusOne = Long.fromNumber(nside-1.); 139 var lOne = Long.fromInt(1); 140 if ( jp.greaterThanOrEqual(lNside) ) 141 jp = nSideMinusOne; 142 if ( jm.greaterThanOrEqual(lNside) ) 143 jm = nSideMinusOne; 144 145 if (z>=0) 146 { 147 return xyf2nest( lNside.subtract(jm).subtract(lOne).toInt(), lNside.subtract(jp).subtract(lOne).toInt(), ntt, order ); 148 } 149 else 150 { 151 return xyf2nest( jp.toInt(), jm.toInt(), ntt+8, order ); 152 } 153 } 154 } 155 156 /**************************************************************************************************************/ 157 158 var HEALPixBase = { 159 compress_bits: function(v){ 160 // raw = v & 0x5555555555555 in place of raw = v & 0x5555555555555555 161 // --> still not resolved, dunno why 162 // 163 164 // in Java implementation mask == 0x5555555555555555 165 // var raw = v & 0x5555555555555; // v & 101010101010101010101010101010101010101010101010101010101010101 166 // // raw>>>15 = 0101010101010101010101010101010101010101010101010 167 // var dec = raw>>>15; 168 // raw |= dec; // 101010101010101111111111111111111111111111111111111111111111111 169 // var raw1 = (raw&0xffff); 170 // var dec2 = raw>>>31; 171 // var raw2 = (dec2&0xffff); 172 173 var longV = Long.fromNumber(v); 174 var longMask = Long.fromNumber(0x5555555555555); 175 var raw = longV.and(longMask); 176 var dec = raw.shiftRightUnsigned(15); 177 raw = raw.or(dec); 178 var raw1 = (raw.and(Long.fromNumber(0xffff))).toInt(); 179 var dec2 = raw.shiftRightUnsigned(32); 180 var raw2 = (dec2.and(Long.fromNumber(0xffff))).toInt(); 181 182 return HealPixTables.ctab[raw1&0xff] | (HealPixTables.ctab[raw1>>>8]<< 4) 183 | (HealPixTables.ctab[raw2&0xff]<<16) | (HealPixTables.ctab[raw2>>>8]<<20); 184 }, 185 186 /** 187 * Function describing a location on the sphere 188 */ 189 fxyf: function(_fx,_fy,_face){ 190 var jr = HealPixTables.jrll[_face] - _fx - _fy; 191 var z = 0; 192 var phi = 0; 193 var sth = 0; 194 var have_sth = false; 195 196 var nr; 197 if (jr<1){ 198 nr = jr; 199 var tmp = nr*nr/3.; 200 z = 1 - tmp; 201 if (z>0.99) { sth=Math.sqrt(tmp*(2.-tmp)); have_sth=true; } 202 } else if (jr>3){ 203 nr = 4-jr; 204 var tmp = nr*nr/3.; 205 z = tmp - 1; 206 if (z<-0.99) { 207 sth=Math.sqrt(tmp*(2.-tmp)); 208 have_sth=true; 209 } 210 } else { 211 nr = 1; 212 z = (2-jr)*2./3.; 213 } 214 215 var tmp=HealPixTables.jpll[_face]*nr+_fx-_fy; 216 if (tmp<0) tmp+=8; 217 if (tmp>=8) tmp-=8; 218 219 phi = (nr<1e-15) ? 0 : (0.5*HALF_PI*tmp)/nr; 220 221 var st = (have_sth) ? sth : Math.sqrt((1.0-z)*(1.0+z)); 222 return [st*Math.cos(phi), st*Math.sin(phi), z]; 223 }, 224 225 /** 226 * Static function 227 * Convert nside to order 228 * (ilog2(nside)) 229 */ 230 nside2order: function(arg){ 231 var res=0; 232 while (arg > 0x0000FFFF) { res+=16; arg>>>=16; } 233 if (arg > 0x000000FF) { res|=8; arg>>>=8; } 234 if (arg > 0x0000000F) { res|=4; arg>>>=4; } 235 if (arg > 0x00000003) { res|=2; arg>>>=2; } 236 if (arg > 0x00000001) { res|=1; } 237 return res; 238 }, 239 240 /** 241 * Returns pixel index of point on sphere 242 * 243 * @param order Tile order 244 * @param lon Longitude 245 * @param lat Latitude 246 */ 247 lonLat2pix: function(order, lon, lat){ 248 var loc = lonLat2ang( lon, lat ); 249 return loc2pix( order, loc[0], loc[1] ); 250 } 251 }; 252 253 /**************************************************************************************************************/ 254 255 return HEALPixBase; 256 257 });