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 });