Changeset 1385 for trunk/psLib/src/astro/psCoord.c
- Timestamp:
- Aug 4, 2004, 1:37:39 PM (22 years ago)
- File:
-
- 1 edited
-
trunk/psLib/src/astro/psCoord.c (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/astro/psCoord.c
r1374 r1385 10 10 * @author George Gusciora, MHPCC 11 11 * 12 * @version $Revision: 1. 8$ $Name: not supported by cvs2svn $13 * @date $Date: 2004-08-04 00:55:17$12 * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2004-08-04 23:37:39 $ 14 14 * 15 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 33 33 { 34 34 if ( out == NULL ) { 35 36 35 out = ( psPlane * ) psAlloc( sizeof( psPlane ) ); 36 } 37 37 out->x = transform->x->coeff[ 0 ][ 0 ] + 38 38 ( transform->x->coeff[ 1 ][ 0 ] * coords->x ) + 39 39 ( transform->x->coeff[ 0 ][ 1 ] * coords->y ); 40 40 41 41 out->y = transform->y->coeff[ 0 ][ 0 ] + 42 42 ( transform->y->coeff[ 1 ][ 0 ] * coords->x ) + 43 43 ( transform->y->coeff[ 0 ][ 1 ] * coords->y ); 44 44 45 45 return ( out ); 46 46 } … … 56 56 { 57 57 if ( out == NULL ) { 58 59 60 58 out = ( psPlane * ) psAlloc( sizeof( psPlane ) ); 59 } 60 61 61 out->x = transform->x->coeff[ 0 ][ 0 ][ 0 ][ 0 ] + 62 62 ( transform->x->coeff[ 1 ][ 0 ][ 0 ][ 0 ] * coords->x ) + … … 64 64 ( transform->x->coeff[ 0 ][ 0 ][ 1 ][ 0 ] * term3 ) + 65 65 ( transform->x->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 ); 66 66 67 67 out->y = transform->y->coeff[ 0 ][ 0 ][ 0 ][ 0 ] + 68 68 ( transform->y->coeff[ 1 ][ 0 ][ 0 ][ 0 ] * coords->x ) + … … 70 70 ( transform->y->coeff[ 0 ][ 0 ][ 1 ][ 0 ] * term3 ) + 71 71 ( transform->y->coeff[ 0 ][ 0 ][ 0 ][ 1 ] * term4 ); 72 72 73 73 return ( out ); 74 74 } … … 81 81 { 82 82 psSphereTransform * tmp = ( psSphereTransform * ) psAlloc( sizeof( psSphereTransform ) ); 83 83 84 84 tmp->sinPhi = sin( NPlat ); 85 85 tmp->cosPhi = cos( NPlat ); 86 86 tmp->Xo = Xo; 87 87 tmp->xo = xo; 88 88 89 89 return ( tmp ); 90 90 } … … 114 114 double y = 0.0; 115 115 double dx = 0.0; 116 116 117 117 if ( out == NULL ) { 118 119 120 118 out = ( psSphere * ) psAlloc( sizeof( psSphere ) ); 119 } 120 121 121 x = coord->r; 122 122 y = coord->d; … … 127 127 cos( y ); 128 128 cosX = cos( y ) * cos( dx ) / cos( y ); 129 129 130 130 out->r = atan2( sinX, cosX ) + transform->Xo; 131 131 out->d = atan2( sinY, cosY ); 132 132 133 133 return ( out ); 134 134 } … … 142 142 double Xo = 0.0; 143 143 double xo = 0.0; 144 144 145 145 return ( psSphereTransformAlloc( phi, Xo, xo ) ); 146 146 } … … 154 154 double Xo = 0.0; 155 155 double xo = 0.0; 156 156 157 157 return ( psSphereTransformAlloc( phi, Xo, xo ) ); 158 158 } … … 178 178 { 179 179 if ( x > 0 ) { 180 return ( atan( y / x ) ); 181 } else if ( ( x == 0 ) && ( y == 0 ) ) { 180 return ( atan( y / x ) ); 181 } else 182 if ( ( x == 0 ) && ( y == 0 ) ) { 182 183 return ( 0.5 * M_PI ); 183 } else if ( ( x == 0 ) && ( y == 0 ) ) { 184 return ( -0.5 * M_PI ); 185 } else if ( ( x == 0 ) && ( y == 0 ) ) { 186 return ( M_PI + atan( y / x ) ); 187 } else if ( ( x == 0 ) && ( y == 0 ) ) { 188 return ( -M_PI + atan( y / x ) ); 189 } 190 184 } else 185 if ( ( x == 0 ) && ( y == 0 ) ) { 186 return ( -0.5 * M_PI ); 187 } else 188 if ( ( x == 0 ) && ( y == 0 ) ) { 189 return ( M_PI + atan( y / x ) ); 190 } else 191 if ( ( x == 0 ) && ( y == 0 ) ) { 192 return ( -M_PI + atan( y / x ) ); 193 } 194 191 195 psAbort( __func__, "Unacceptable range for (arg(%f, %f).\n", x, y ); 192 196 return ( 0.0 ); … … 201 205 float alpha = 0.0; 202 206 psPlane *tmp = ( psPlane * ) psAlloc( sizeof( psPlane ) ); 203 207 204 208 if ( projection->type == PS_PROJ_TAN ) { 205 R = cot( coord->r ) * ( 180.0 / M_PI ); 206 tmp->x = R * sin( coord->d ); 207 tmp->y = R * cos( coord->d ); 208 209 } else if ( projection->type == PS_PROJ_SIN ) { 209 R = cot( coord->r ) * ( 180.0 / M_PI ); 210 tmp->x = R * sin( coord->d ); 211 tmp->y = R * cos( coord->d ); 212 213 } else 214 if ( projection->type == PS_PROJ_SIN ) { 210 215 R = cos( coord->r ) * ( 180.0 / M_PI ); 211 216 tmp->x = R * sin( coord->d ); 212 217 tmp->y = R * cos( coord->d ); 213 214 } else if ( projection->type == PS_PROJ_CAR ) { 215 tmp->x = coord->d; 216 tmp->y = coord->r; 217 218 } else if ( projection->type == PS_PROJ_MER ) { 219 tmp->x = coord->d; 220 tmp->y = log( tan( 45.0 + ( 0.5 * coord->r ) ) ) * 180.0 / M_PI; 221 222 } else if ( projection->type == PS_PROJ_AIT ) { 223 alpha = 1.0 / ( ( 180.0 / M_PI ) * 224 sqrt( 1.0 + ( cos( coord->r ) * cos( 0.5 * coord->d ) * 0.5 ) ) ); 225 226 tmp->x = 2.0 * alpha * cos( coord->r ) * sin( 0.5 * coord->d ); 227 tmp->y = alpha * sin( coord->d ); 228 229 } else if ( projection->type == PS_PROJ_PAR ) { 230 psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" ); 231 232 } else if ( projection->type == PS_PROJ_GLS ) { 233 psAbort( __func__, "The projection type PS_PROJ_GLS is undefined.\n" ); 234 } 235 218 219 } else 220 if ( projection->type == PS_PROJ_CAR ) { 221 tmp->x = coord->d; 222 tmp->y = coord->r; 223 224 } else 225 if ( projection->type == PS_PROJ_MER ) { 226 tmp->x = coord->d; 227 tmp->y = log( tan( 45.0 + ( 0.5 * coord->r ) ) ) * 180.0 / M_PI; 228 229 } else 230 if ( projection->type == PS_PROJ_AIT ) { 231 alpha = 1.0 / ( ( 180.0 / M_PI ) * 232 sqrt( 1.0 + ( cos( coord->r ) * cos( 0.5 * coord->d ) * 0.5 ) ) ); 233 234 tmp->x = 2.0 * alpha * cos( coord->r ) * sin( 0.5 * coord->d ); 235 tmp->y = alpha * sin( coord->d ); 236 237 } else 238 if ( projection->type == PS_PROJ_PAR ) { 239 psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" ); 240 241 } else 242 if ( projection->type == PS_PROJ_GLS ) { 243 psAbort( __func__, "The projection type PS_PROJ_GLS is undefined.\n" ); 244 } 245 236 246 return ( tmp ); 237 247 } … … 247 257 float chu2 = 0.0; 248 258 psSphere *tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) ); 249 259 250 260 if ( projection->type == PS_PROJ_TAN ) { 251 R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) ); 252 tmp->d = arg( -coord->y, coord->x ); 253 tmp->r = atan( 180.0 / ( R * M_PI ) ); 254 255 } else if ( projection->type == PS_PROJ_SIN ) { 261 R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) ); 262 tmp->d = arg( -coord->y, coord->x ); 263 tmp->r = atan( 180.0 / ( R * M_PI ) ); 264 265 } else 266 if ( projection->type == PS_PROJ_SIN ) { 256 267 R = sqrt( ( coord->x * coord->x ) + ( coord->y * coord->y ) ); 257 268 tmp->d = arg( -coord->y, coord->x ); 258 269 tmp->r = acos( ( R * M_PI ) / 180.0 ); 259 260 } else if ( projection->type == PS_PROJ_CAR ) { 261 tmp->d = coord->x; 262 tmp->r = coord->y; 263 264 } else if ( projection->type == PS_PROJ_MER ) { 265 tmp->d = coord->x; 266 tmp->r = ( 2.0 * atan( exp( ( coord->y * M_PI / 180.0 ) ) ) ) - 180.0; 267 268 } else if ( projection->type == PS_PROJ_AIT ) { 269 chu1 = ( coord->x * M_PI ) / 720.0; 270 chu1 *= chu1; 271 chu2 = ( coord->y * M_PI ) / 360.0; 272 chu2 *= chu2; 273 chu = sqrt( 1.0 - chu1 - chu2 ); 274 tmp->d = 2.0 * arg( ( 2.0 * chu * chu ) - 1.0, 275 ( coord->x * chu * M_PI ) / 360.0 ); 276 tmp->r = asin( ( coord->y * chu * M_PI ) / 180.0 ); 277 278 } else if ( projection->type == PS_PROJ_PAR ) { 279 psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" ); 280 281 } else if ( projection->type == PS_PROJ_GLS ) { 282 psAbort( __func__, "The projection type PS_PROJ_GLG is undefined.\n" ); 283 } 284 270 271 } else 272 if ( projection->type == PS_PROJ_CAR ) { 273 tmp->d = coord->x; 274 tmp->r = coord->y; 275 276 } else 277 if ( projection->type == PS_PROJ_MER ) { 278 tmp->d = coord->x; 279 tmp->r = ( 2.0 * atan( exp( ( coord->y * M_PI / 180.0 ) ) ) ) - 180.0; 280 281 } else 282 if ( projection->type == PS_PROJ_AIT ) { 283 chu1 = ( coord->x * M_PI ) / 720.0; 284 chu1 *= chu1; 285 chu2 = ( coord->y * M_PI ) / 360.0; 286 chu2 *= chu2; 287 chu = sqrt( 1.0 - chu1 - chu2 ); 288 tmp->d = 2.0 * arg( ( 2.0 * chu * chu ) - 1.0, 289 ( coord->x * chu * M_PI ) / 360.0 ); 290 tmp->r = asin( ( coord->y * chu * M_PI ) / 180.0 ); 291 292 } else 293 if ( projection->type == PS_PROJ_PAR ) { 294 psAbort( __func__, "The projection type PS_PROJ_PAR is undefined.\n" ); 295 296 } else 297 if ( projection->type == PS_PROJ_GLS ) { 298 psAbort( __func__, "The projection type PS_PROJ_GLG is undefined.\n" ); 299 } 300 285 301 return ( tmp ); 286 302 } … … 298 314 double tmpR = 0.0; 299 315 double tmpD = 0.0; 300 316 301 317 if ( mode == PS_LINEAR ) { 302 // XXX: I have no idea how to construct this. Maybe project both 303 // sperical positions onto the plane, set the origin at one of the 304 // points on the plane, then deproject? 305 306 // XXX: Do I need to somehow scale this projection? 307 // project position1? Will it project to (0.0, 0.0)? 308 proj.R = position1->r; 309 proj.D = position1->d; 310 proj.Xs = 1.0; 311 proj.Ys = 1.0; 312 proj.type = PS_PROJ_TAN; 313 314 // lin = psProject(position2, proj); 315 // tmp = psDeproject(lin, proj); 316 317 // XXX: Do we need to convert units in tmp? 318 return ( tmp ); 319 } else if ( mode == PS_SPHERICAL ) { 318 // XXX: I have no idea how to construct this. Maybe project both 319 // sperical positions onto the plane, set the origin at one of the 320 // points on the plane, then deproject? 321 322 // XXX: Do I need to somehow scale this projection? 323 // project position1? Will it project to (0.0, 0.0)? 324 proj.R = position1->r; 325 proj.D = position1->d; 326 proj.Xs = 1.0; 327 proj.Ys = 1.0; 328 proj.type = PS_PROJ_TAN; 329 330 // lin = psProject(position2, proj); 331 // tmp = psDeproject(lin, proj); 332 333 // XXX: Do we need to convert units in tmp? 334 return ( tmp ); 335 } else 336 if ( mode == PS_SPHERICAL ) { 320 337 tmpR = position2->r - position1->r; 321 338 tmpD = position2->d - position1->d; 322 339 323 340 if ( unit == PS_ARCSEC ) { 324 tmpR = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI; 325 tmpD = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI; 326 } else if ( unit == PS_ARCMIN ) { 341 tmpR = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI; 342 tmpD = ( tmpR * 180.0 * 60.0 * 60.0 ) / M_PI; 343 } else 344 if ( unit == PS_ARCMIN ) { 327 345 tmpR = ( tmpR * 180.0 * 60.0 ) / M_PI; 328 346 tmpD = ( tmpR * 180.0 * 60.0 ) / M_PI; 329 } else if ( unit == PS_DEGREE ) { 330 tmpR = ( tmpR * 180.0 ) / M_PI; 331 tmpD = ( tmpR * 180.0 ) / M_PI; 332 } else if ( unit == PS_RADIAN ) {} 333 else { 334 psAbort( __func__, "Unknown offset unit: 0x%x\n", unit ); 335 } 336 347 } else 348 if ( unit == PS_DEGREE ) { 349 tmpR = ( tmpR * 180.0 ) / M_PI; 350 tmpD = ( tmpR * 180.0 ) / M_PI; 351 } else 352 if ( unit == PS_RADIAN ) {} 353 else { 354 psAbort( __func__, "Unknown offset unit: 0x%x\n", unit ); 355 } 356 337 357 tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) ); 338 358 tmp->r = tmpR; … … 361 381 double tmpR = 0.0; 362 382 double tmpD = 0.0; 363 383 364 384 if ( mode == PS_LINEAR ) { 365 proj.R = position->r; 366 proj.D = position->d; 367 proj.Xs = 1.0; 368 proj.Ys = 1.0; 369 proj.type = PS_PROJ_TAN; 370 371 lin.x = offset->r; 372 lin.y = offset->d; 373 374 tmp = psDeproject( &lin, &proj ); 375 return ( tmp ); 376 377 } else if ( mode == PS_SPHERICAL ) { 385 proj.R = position->r; 386 proj.D = position->d; 387 proj.Xs = 1.0; 388 proj.Ys = 1.0; 389 proj.type = PS_PROJ_TAN; 390 391 lin.x = offset->r; 392 lin.y = offset->d; 393 394 tmp = psDeproject( &lin, &proj ); 395 return ( tmp ); 396 397 } else 398 if ( mode == PS_SPHERICAL ) { 378 399 if ( unit == PS_ARCSEC ) { 379 tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 * 60.0 ); 380 tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 * 60.0 ); 381 } else if ( unit == PS_ARCMIN ) { 400 tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 * 60.0 ); 401 tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 * 60.0 ); 402 } else 403 if ( unit == PS_ARCMIN ) { 382 404 tmpR = ( M_PI * offset->r ) / ( 180.0 * 60.0 ); 383 405 tmpD = ( M_PI * offset->d ) / ( 180.0 * 60.0 ); 384 } else if ( unit == PS_DEGREE ) { 385 tmpR = ( M_PI * offset->r ) / ( 180.0 ); 386 tmpD = ( M_PI * offset->d ) / ( 180.0 ); 387 } else if ( unit == PS_RADIAN ) { 388 tmpR = offset->r; 389 tmpD = offset->d; 390 } else { 391 psAbort( __func__, "Unknown offset unit: 0x%x\n", unit ); 392 } 393 406 } else 407 if ( unit == PS_DEGREE ) { 408 tmpR = ( M_PI * offset->r ) / ( 180.0 ); 409 tmpD = ( M_PI * offset->d ) / ( 180.0 ); 410 } else 411 if ( unit == PS_RADIAN ) { 412 tmpR = offset->r; 413 tmpD = offset->d; 414 } else { 415 psAbort( __func__, "Unknown offset unit: 0x%x\n", unit ); 416 } 417 394 418 tmp = ( psSphere * ) psAlloc( sizeof( psSphere ) ); 395 419 tmp->r = position->r + tmpR; … … 397 421 tmp->rErr = 0.0; 398 422 tmp->dErr = 0.0; 399 423 400 424 // XXX: wrap tmp->r and tmp->d to the allowed range (-PI to PI) 401 425 // and (0 to 2*PI).
Note:
See TracChangeset
for help on using the changeset viewer.
