001 /*
002 * This file is part of McIDAS-V
003 *
004 * Copyright 2007-2013
005 * Space Science and Engineering Center (SSEC)
006 * University of Wisconsin - Madison
007 * 1225 W. Dayton Street, Madison, WI 53706, USA
008 * https://www.ssec.wisc.edu/mcidas
009 *
010 * All Rights Reserved
011 *
012 * McIDAS-V is built on Unidata's IDV and SSEC's VisAD libraries, and
013 * some McIDAS-V source code is based on IDV and VisAD source code.
014 *
015 * McIDAS-V is free software; you can redistribute it and/or modify
016 * it under the terms of the GNU Lesser Public License as published by
017 * the Free Software Foundation; either version 3 of the License, or
018 * (at your option) any later version.
019 *
020 * McIDAS-V is distributed in the hope that it will be useful,
021 * but WITHOUT ANY WARRANTY; without even the implied warranty of
022 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
023 * GNU Lesser Public License for more details.
024 *
025 * You should have received a copy of the GNU Lesser Public License
026 * along with this program. If not, see http://www.gnu.org/licenses.
027 */
028
029 package edu.wisc.ssec.mcidasv.data.adde.sgp4;
030
031 //import name.gano.astro.AstroConst;
032 //import name.gano.astro.MathUtils;
033
034 /**
035 *
036 * @author Shawn
037 */
038 public class J2kCoordinateConversion
039 {
040 public static enum Opt
041 {
042 e80,
043 e96,
044 e00a,
045 e00b
046 }
047
048 public static enum Direction
049 {
050 to,
051 from
052 }
053
054 /* -----------------------------------------------------------------------------
055 *
056 * function teme_j2k
057 *
058 * this function transforms a vector between the true equator mean equinox system,
059 * (teme) and the mean equator mean equinox (j2000) system.
060 *
061 * author : david vallado 719-573-2600 25 jun 2002
062 *
063 * revisions
064 * vallado - add order 29 sep 2002
065 * vallado - conversion to c++ 21 feb 2005
066 *
067 * inputs description range / units
068 * rteme - position vector of date
069 * true equator, mean equinox km
070 * vteme - velocity vector of date
071 * true equator, mean equinox km/s
072 * ateme - acceleration vector of date
073 * true equator, mean equinox km/s2
074 * direct - direction of transfer eFrom, 'TOO '
075 * iau80rec - record containing the iau80 constants rad
076 * ttt - julian centuries of tt centuries
077 * eqeterms - number of terms for eqe 0, 2
078 * opt1 - option for processing a - complete nutation
079 * b - truncated nutation
080 * c - truncated transf matrix
081 *
082 * outputs :
083 * rj2k - position vector j2k km
084 * vj2k - velocity vector j2k km/s
085 * aj2k - acceleration vector j2k km/s2
086 *
087 * locals :
088 * prec - matrix for mod - j2k
089 * nutteme - matrix for mod - teme - an approximation for nutation
090 * tm - combined matrix for teme
091 *
092 * coupling :
093 * precess - rotation for precession j2k - mod
094 * truemean - rotation for truemean j2k - teme
095 *
096 * references :
097 * vallado 2007, 236
098 * --------------------------------------------------------------------------- */
099 //returns the transformation matrix
100 public static double[][] teme_j2k
101 (
102 //double rteme[3], double vteme[3], double ateme[3],
103 Direction direct,
104 //double rj2k[3], double vj2k[3], double aj2k[3],
105 //iau80data& iau80rec,
106 double ttt, int order, int eqeterms, char optteme
107 )
108 {
109 //std::vector< std::vector<double> > prec, nutteme, temp[3][3], tempmat, nuttemep, precp;
110 double[][] tempmat = new double[3][3];
111
112 // double psia, wa, epsa, chia, deltapsi, trueeps, meaneps,
113 // omega, thetasa;
114 // double[] omegaearth = new double[3];
115 // double[] omgxr = new double[3];
116 // double[] omgxomgxr = new double[3];
117 // double[] omgxv = new double[3];
118 // double[] tempvec1 = new double[3];
119 // double[] tempvec = new double[3];
120
121 double[][] prec = precess ( ttt, Opt.e80);//, psia,wa,epsa,chia );
122 double[][] nutteme = truemean( ttt,order,eqeterms,optteme );
123
124 if (direct == Direction.to)
125 {
126 tempmat = matmult( prec, nutteme, 3, 3, 3 );
127 // double[] rj2k = matvecmult( tempmat, rteme);
128 // double[] vj2k = matvecmult( tempmat, vteme);
129 // double[] aj2k = matvecmult( tempmat, ateme);
130 }
131 else
132 {
133 double[][] nuttemep = mattrans(nutteme, 3, 3 );
134 double[][] precp = mattrans(prec, 3, 3 );
135
136 tempmat = matmult( nuttemep, precp, 3, 3, 3 );
137
138 // double[] rteme = matvecmult(tempmat, rj2k );
139 // double[] vteme = matvecmult(tempmat, vj2k);
140 // double[] ateme = matvecmult(tempmat, aj2k);
141 }
142
143 return tempmat;
144 } // procedure teme_j2k
145
146 /* -----------------------------------------------------------------------------
147 *
148 * function precess
149 *
150 * this function calulates the transformation matrix that accounts for the effects
151 * of precession. both the 1980 and 2000 theories are handled. note that the
152 * required parameters differ a little.
153 *
154 * author : david vallado 719-573-2600 25 jun 2002
155 *
156 * revisions
157 * vallado - conversion to c++ 21 feb 2005
158 * vallado - misc updates, nomenclature, etc 23 nov 2005
159 *
160 * inputs description range / units
161 * ttt - julian centuries of tt
162 * opt - method option e00a, e00b, e96, e80
163 *
164 * outputs :
165 * prec - transformation matrix for mod - j2000 (80 only)
166 * psia - cannonical precession angle rad (00 only)
167 * wa - cannonical precession angle rad (00 only)
168 * epsa - cannonical precession angle rad (00 only)
169 * chia - cannonical precession angle rad (00 only)
170 * prec - matrix converting from "mod" to gcrf
171 *
172 * locals :
173 * zeta - precession angle rad
174 * z - precession angle rad
175 * theta - precession angle rad
176 * oblo - obliquity value at j2000 epoch "//
177 *
178 * coupling :
179 * none -
180 *
181 * references :
182 * vallado 2007, 217, 228
183 * --------------------------------------------------------------------------- */
184 public static double[][] precess
185 (
186 double ttt, Opt opt//,
187 //double& psia, double& wa, double& epsa, double& chia,
188 //std::vector< std::vector<double> > &prec
189 )
190 {
191 double[][] prec = new double[3][3];//.resize(3); // rows
192 //for (std::vector< std::vector<double> >::iterator it=prec.begin(); it != prec.end();++it)
193 // it->resize(3);
194 //std::vector< std::vector<double> > p1, p2, p3, p4, tr1, tr2;
195 double[][] p1= new double[3][3];
196 double[][] p2= new double[3][3];
197 double[][] p3= new double[3][3];
198 double[][] p4= new double[3][3];
199 double[][] tr1= new double[3][3];
200 double[][] tr2= new double[3][3];
201
202 // since not returning these for now
203 double psia, wa, epsa, chia;
204
205 double convrt, zeta, theta, z, coszeta, sinzeta, costheta, sintheta,
206 cosz, sinz, oblo;
207 //char iauhelp;
208 //sethelp(iauhelp, ' ');
209
210 convrt = Math.PI / (180.0 * 3600.0);
211
212 // ------------------- iau 77 precession angles --------------------
213 if ((opt == Opt.e80) | (opt == Opt.e96))
214 {
215 oblo = 84381.448; // "
216 psia = (( - 0.001147 * ttt - 1.07259 ) * ttt + 5038.7784 ) * ttt; // "
217 wa = (( - 0.007726 * ttt + 0.05127 ) * ttt ) + oblo;
218 epsa = (( 0.001813 * ttt - 0.00059 ) * ttt - 46.8150 ) * ttt + oblo;
219 chia = (( - 0.001125 * ttt - 2.38064 ) * ttt + 10.5526 ) * ttt;
220
221 zeta = (( 0.017998 * ttt + 0.30188 ) * ttt + 2306.2181 ) * ttt; // "
222 theta= (( - 0.041833 * ttt - 0.42665 ) * ttt + 2004.3109 ) * ttt;
223 z = (( 0.018203 * ttt + 1.09468 ) * ttt + 2306.2181 ) * ttt;
224 }
225 // ------------------ iau 00 precession angles -------------------
226 else
227 {
228 oblo = 84381.406; // "
229 psia = (((( -0.0000000951 * ttt + 0.000132851 ) * ttt - 0.00114045 ) * ttt - 1.0790069 ) * ttt + 5038.481507 ) * ttt; // "
230 wa = (((( 0.0000003337 * ttt - 0.000000467 ) * ttt - 0.00772503 ) * ttt + 0.0512623 ) * ttt - 0.025754 ) * ttt + oblo;
231 epsa = (((( -0.0000000434 * ttt - 0.000000576 ) * ttt + 0.00200340 ) * ttt - 0.0001831 ) * ttt - 46.836769 ) * ttt + oblo;
232 chia = (((( - 0.0000000560 * ttt + 0.000170663 ) * ttt - 0.00121197 ) * ttt - 2.3814292 ) * ttt + 10.556403 ) * ttt;
233
234 zeta = (((( - 0.0000003173 * ttt - 0.000005971 ) * ttt + 0.01801828 ) * ttt + 0.2988499 ) * ttt + 2306.083227 ) * ttt + 2.650545; // "
235 theta= (((( - 0.0000001274 * ttt - 0.000007089 ) * ttt - 0.04182264 ) * ttt - 0.4294934 ) * ttt + 2004.191903 ) * ttt;
236 z = (((( 0.0000002904 * ttt - 0.000028596 ) * ttt + 0.01826837 ) * ttt + 1.0927348 ) * ttt + 2306.077181 ) * ttt - 2.650545;
237 }
238
239 // convert units to rad
240 psia = psia * convrt;
241 wa = wa * convrt;
242 oblo = oblo * convrt;
243 epsa = epsa * convrt;
244 chia = chia * convrt;
245
246 zeta = zeta * convrt;
247 theta= theta * convrt;
248 z = z * convrt;
249
250 if ((opt == Opt.e80) | (opt == Opt.e96))
251 {
252 coszeta = Math.cos(zeta);
253 sinzeta = Math.sin(zeta);
254 costheta = Math.cos(theta);
255 sintheta = Math.sin(theta);
256 cosz = Math.cos(z);
257 sinz = Math.sin(z);
258
259 // ----------------- form matrix mod to gcrf ------------------
260 prec[0][0] = coszeta * costheta * cosz - sinzeta * sinz;
261 prec[0][1] = coszeta * costheta * sinz + sinzeta * cosz;
262 prec[0][2] = coszeta * sintheta;
263 prec[1][0] = -sinzeta * costheta * cosz - coszeta * sinz;
264 prec[1][1] = -sinzeta * costheta * sinz + coszeta * cosz;
265 prec[1][2] = -sinzeta * sintheta;
266 prec[2][0] = -sintheta * cosz;
267 prec[2][1] = -sintheta * sinz;
268 prec[2][2] = costheta;
269
270 // ----------------- do rotations instead ----------------------
271 // p1 = rot3mat( z );
272 // p2 = rot2mat( -theta );
273 // p3 = rot3mat( zeta );
274 // prec = p3 * p2*p1;
275 }
276 else
277 {
278 p1 = rot3mat( -chia );
279 p2 = rot1mat( wa );
280 p3 = rot3mat( psia );
281 p4 = rot1mat( -oblo );
282
283 tr1 = matmult( p4 , p3, 3, 3, 3 );
284 tr2 = matmult( tr1, p2, 3, 3, 3 );
285 prec = matmult( tr2, p1, 3, 3, 3 );
286
287 }
288
289 // if (iauhelp == 'y')
290 // {
291 // printf("pr %11.7f %11.7f %11.7f %11.7fdeg \n",psia * 180/pi,wa * 180/pi,epsa * 180/pi,chia * 180/pi );
292 // printf("pr %11.7f %11.7f %11.7fdeg \n",zeta * 180/pi,theta * 180/pi,z * 180/pi );
293 // }
294
295 return prec;
296
297 } // procedure precess
298
299
300 /* -----------------------------------------------------------------------------
301 *
302 * procedure mattrans
303 *
304 * this procedure finds the transpose of a matrix.
305 *
306 * author : david vallado 719-573-2600 1 mar 2001
307 *
308 * inputs description range / units
309 * mat1 - matrix number 1
310 * mat1r - matrix number 1 rows
311 * mat1c - matrix number 1 columns
312 *
313 * outputs :
314 * mat2 - matrix result of transpose mat2
315 *
316 * locals :
317 * row - row index
318 * col - column index
319 *
320 * coupling :
321 *
322 * --------------------------------------------------------------------------- */
323
324 public static double[][] mattrans
325 (
326 double[][] mat1,
327 //double[][] mat2,
328 // double mat1[3][3],
329 // double mat2[3][3],
330 int mat1r, int mat1c
331 )
332 {
333 int row,col;
334
335 double[][] mat2 = new double[mat1c][mat1r];
336 //mat2.resize(mat1c); // rows
337 //for (std::vector< std::vector<double> >::iterator it=mat2.begin(); it != mat2.end();++it)
338 // it->resize(mat1r);
339
340 for (row = 0; row < mat1r; row++)
341 {
342 for (col = 0; col < mat1c; col++)
343 mat2[col][row] = mat1[row][col];
344 }
345
346 return mat2;
347 }
348
349
350 /* -----------------------------------------------------------------------------
351 *
352 * procedure matmult
353 *
354 * this procedure multiplies two matricies up to 10x10 together.
355 *
356 * author : david vallado 719-573-2600 7 dec 2007
357 *
358 * inputs description range / units
359 * mat1 - matrix number 1
360 * mat2 - matrix number 2
361 * mat1r - matrix number 1 rows
362 * mat1c - matrix number 1 columns
363 * mat2c - matrix number 2 columns
364 *
365 * outputs :
366 * mat3 - matrix result of mat1 * mat2 of size mat1r x mat2c
367 *
368 * locals :
369 * row - row index
370 * col - column index
371 * ktr - index
372 *
373 * coupling :
374 *
375 * --------------------------------------------------------------------------- */
376
377 public static double[][] matmult
378 (
379 double[][] mat1,
380 double[][] mat2,
381 // std::vector< std::vector<double> > &mat3,
382 // double mat1[3][3],
383 // double mat2[3][3],
384 // double mat3[3][3],
385 int mat1r, int mat1c, int mat2c
386 )
387 {
388 double[][] mat3 = new double[mat1r][mat2c];
389 int row,col,ktr;
390 // specify the actual sizes
391 // mat3.resize(mat1r); // rows
392 // for (std::vector< std::vector<double> >::iterator it=mat3.begin(); it != mat3.end();++it)
393 // it->resize(mat2c);
394
395
396 for (row = 0; row < mat1r; row++)
397 {
398 for (col = 0; col < mat2c; col++)
399 {
400 mat3[row][col]= 0.0;
401 for (ktr = 0; ktr < mat1c; ktr ++)
402 mat3[row][col]= mat3[row][col] + mat1[row][ktr] * mat2[ktr][col];
403 }
404 }
405
406 return mat3;
407 } //matmult
408
409 /* -----------------------------------------------------------------------------
410 *
411 * procedure matvecmult
412 *
413 * this procedure multiplies a 3x3 matrix and a 3x1 vector together.
414 *
415 * author : david vallado 719-573-2600 1 mar 2001
416 *
417 * inputs description range / units
418 * mat - 3 x 3 matrix
419 * vec - vector
420 *
421 * outputs :
422 * vecout - vector result of mat * vec
423 *
424 * locals :
425 * row - row index
426 * col - column index
427 * ktr - index
428 *
429 * coupling :
430 *
431 * --------------------------------------------------------------------------- */
432
433 public static double[] matvecmult
434 (
435 double[][] mat,
436 // double mat[3][3],
437 double[] vec//,
438 //double[] vecout
439 )
440 {
441 int row,col,ktr;
442 double[] vecout = new double[3];
443
444 for (row = 0; row <= 2; row++)
445 {
446 vecout[row]= 0.0;
447 for (ktr = 0; ktr <= 2; ktr++)
448 vecout[row]= vecout[row] + mat[row][ktr] * vec[ktr];
449 }
450
451 return vecout;
452 }
453
454 /* -----------------------------------------------------------------------------
455 *
456 * rotimat
457 *
458 * this function sets up a rotation matrix for an input angle about the first
459 * axis.
460 *
461 * author : david vallado 719-573-2600 10 jan 2003
462 *
463 * revisions
464 * -
465 *
466 * inputs description range / units
467 * xval - angle of rotation rad
468 *
469 * outputs :
470 * outmat - matrix result
471 *
472 * locals :
473 * c - cosine of the angle xval
474 * s - sine of the angle xval
475 *
476 * coupling :
477 * none.
478 *
479 * --------------------------------------------------------------------------- */
480 public static double[][] rot1mat
481 (
482 double xval//,
483 //std::vector< std::vector<double> > &outmat
484 // double outmat[3][3]
485 )
486 {
487 double[][] outmat = new double[3][3];//outmat.resize(3); // rows
488 //for (std::vector< std::vector<double> >::iterator it=outmat.begin(); it != outmat.end();++it)
489 // it->resize(3);
490 double c, s;
491 c= Math.cos( xval );
492 s= Math.sin( xval );
493
494 outmat[0][0]= 1.0;
495 outmat[0][1]= 0.0;
496 outmat[0][2]= 0.0;
497
498 outmat[1][0]= 0.0;
499 outmat[1][1]= c;
500 outmat[1][2]= s;
501
502 outmat[2][0]= 0.0;
503 outmat[2][1]= -s;
504 outmat[2][2]= c;
505
506 return outmat;
507 }
508
509 public static double[][] rot3mat
510 (
511 double xval//,
512 //std::vector< std::vector<double> > &outmat
513 // double outmat[3][3]
514 )
515 {
516 double[][] outmat = new double[3][3];//.resize(3); // rows
517 //for (std::vector< std::vector<double> >::iterator it=outmat.begin(); it != outmat.end();++it)
518 // it->resize(3);
519 double c, s;
520 c= Math.cos( xval );
521 s= Math.sin( xval );
522
523 outmat[0][0]= c;
524 outmat[0][1]= s;
525 outmat[0][2]= 0.0;
526
527 outmat[1][0]= -s;
528 outmat[1][1]= c;
529 outmat[1][2]= 0.0;
530
531 outmat[2][0]= 0.0;
532 outmat[2][1]= 0.0;
533 outmat[2][2]= 1.0;
534
535 return outmat;
536 }
537
538 /* -----------------------------------------------------------------------------
539 *
540 * function truemean
541 *
542 * this function forms the transformation matrix to go between the
543 * norad true equator mean equinox of date and the mean equator mean equinox
544 * of date (gcrf). the results approximate the effects of nutation and
545 * precession.
546 *
547 * author : david vallado 719-573-2600 25 jun 2002
548 *
549 * revisions
550 * vallado - fixes to order 29 sep 2002
551 * vallado - fixes to all options 6 may 2003
552 * vallado - conversion to c++ 21 feb 2005
553 *
554 * inputs description range / units
555 * ttt - julian centuries of tt
556 * order - number of terms for nutation 4, 50, 106, ...
557 * eqeterms - number of terms for eqe 0, 2
558 * opt1 - option for procesMath.sing a - complete nutation
559 * b - truncated nutation
560 * c - truncated transf matrix
561 * iau80rec - record containing the iau80 constants rad
562 *
563 * outputs :
564 * nutteme - matrix for mod - teme - an approximation for nutation
565 *
566 * locals :
567 * prec - matrix for mod - j2000
568 * tm - combined matrix for teme
569 * l - rad
570 * ll - rad
571 * f - rad
572 * d - rad
573 * omega - rad
574 * deltapsi - nutation angle rad
575 * deltaeps - change in obliquity rad
576 * eps - mean obliquity of the ecliptic rad
577 * trueeps - true obliquity of the ecliptic rad
578 * meaneps - mean obliquity of the ecliptic rad
579 *
580 * coupling :
581 *
582 *
583 * references :
584 * vallado 2007, 236
585 * --------------------------------------------------------------------------- */
586
587 public static double[][] truemean
588 (
589 double ttt, int order, int eqeterms, char opt//,
590 //iau80data& iau80rec,
591 //std::vector< std::vector<double> > &nutteme
592 )
593 {
594 //nutteme.resize(3); // rows
595 //for (std::vector< std::vector<double> >::iterator it=nutteme.begin(); it != nutteme.end();++it)
596 // it->resize(3);
597 double[][] nutteme = new double[3][3];
598
599 double deg2rad, l, l1, f, d, omega,
600 cospsi, sinpsi, coseps, sineps, costrueeps, sintrueeps, meaneps,
601 deltapsi, deltaeps, trueeps;
602 int i;
603 double[][] nut = new double[3][3];
604 double[][] st = new double[3][3];
605 double tempval, jdttt, eqe;
606
607 deg2rad = Math.PI/180.0;
608
609 // ---- determine coefficients for iau 1980 nutation theory ----
610 meaneps = ((0.001813 * ttt - 0.00059 ) * ttt -46.8150 ) * ttt + 84381.448;
611 meaneps = meaneps/3600.0 % 360.0 ;
612 meaneps = meaneps * deg2rad;
613
614 l = (((( 0.064 ) * ttt + 31.310 ) * ttt + 1717915922.6330 ) * ttt )
615 / 3600.0 + 134.96298139;
616 l1 = ((((- 0.012 ) * ttt - 0.577 ) * ttt + 129596581.2240 ) * ttt )
617 / 3600.0 + 357.52772333;
618 f = (((( 0.011 ) * ttt - 13.257 ) * ttt + 1739527263.1370 ) * ttt )
619 / 3600.0 + 93.27191028;
620 d = (((( 0.019 ) * ttt - 6.891 ) * ttt + 1602961601.3280 ) * ttt )
621 / 3600.0 + 297.85036306;
622 omega= (((( 0.008 ) * ttt + 7.455 ) * ttt - 6962890.5390 ) * ttt )
623 / 3600.0 + 125.04452222;
624
625 l = ( l%360.0 ) * deg2rad;
626 l1 = ( l1%360.0 ) * deg2rad;
627 f = ( f%360.0 ) * deg2rad;
628 d = ( d%360.0 ) * deg2rad;
629 omega= ( omega%360.0 ) * deg2rad;
630
631 deltapsi= 0.0;
632 deltaeps= 0.0;
633 for (i= 1; i <= order; i++) // the eqeterms in nut80.dat are already sorted
634 {
635 tempval= iar80(i,1) * l + iar80(i,2) * l1 + iar80(i,3) * f +
636 iar80(i,4) * d + iar80(i,5) * omega;
637 deltapsi= deltapsi + (rar80(i,1)+rar80(i,2) * ttt) * Math.sin( tempval );
638 deltaeps= deltaeps + (rar80(i,3)+rar80(i,4) * ttt) * Math.cos( tempval );
639 }
640
641 // --------------- find nutation parameters --------------------
642 deltapsi = ( deltapsi%360.0 ) * deg2rad;
643 deltaeps = ( deltaeps%360.0 ) * deg2rad;
644 trueeps = meaneps + deltaeps;
645
646 cospsi = Math.cos(deltapsi);
647 sinpsi = Math.sin(deltapsi);
648 coseps = Math.cos(meaneps);
649 sineps = Math.sin(meaneps);
650 costrueeps = Math.cos(trueeps);
651 sintrueeps = Math.sin(trueeps);
652
653 jdttt = ttt * 36525.0 + 2451545.0;
654 // small disconnect with ttt instead of ut1
655 if ((jdttt > 2450449.5 ) && (eqeterms > 0))
656 eqe= deltapsi * Math.cos(meaneps)
657 + 0.00264 * Math.PI /(3600 * 180) * Math.sin(omega)
658 + 0.000063 * Math.PI /(3600 * 180) * Math.sin(2.0 * omega);
659 else
660 eqe= deltapsi * Math.cos(meaneps);
661
662 nut[0][0] = cospsi;
663 nut[0][1] = costrueeps * sinpsi;
664 if (opt == 'b')
665 nut[0][1] = 0.0;
666 nut[0][2] = sintrueeps * sinpsi;
667 nut[1][0] = -coseps * sinpsi;
668 if (opt == 'b')
669 nut[1][0] = 0.0;
670 nut[1][1] = costrueeps * coseps * cospsi + sintrueeps * sineps;
671 nut[1][2] = sintrueeps * coseps * cospsi - sineps * costrueeps;
672 nut[2][0] = -sineps * sinpsi;
673 nut[2][1] = costrueeps * sineps * cospsi - sintrueeps * coseps;
674 nut[2][2] = sintrueeps * sineps * cospsi + costrueeps * coseps;
675
676 st[0][0] = Math.cos(eqe);
677 st[0][1] = -Math.sin(eqe);
678 st[0][2] = 0.0;
679 st[1][0] = Math.sin(eqe);
680 st[1][1] = Math.cos(eqe);
681 st[1][2] = 0.0;
682 st[2][0] = 0.0;
683 st[2][1] = 0.0;
684 st[2][2] = 1.0;
685
686 //matmult( st, nut, nutteme, 3, 3, 3 );
687 nutteme = MathUtils.mult(st, nut);
688
689 if (opt == 'c')
690 {
691 nutteme[0][0] = 1.0;
692 nutteme[0][1] = 0.0;
693 nutteme[0][2] = deltapsi * sineps;
694 nutteme[1][0] = 0.0;
695 nutteme[1][1] = 1.0;
696 nutteme[1][2] = deltaeps;
697 nutteme[2][0] = -deltapsi * sineps;
698 nutteme[2][1] = -deltaeps;
699 nutteme[2][2] = 1.0;
700 }
701
702 // tm = nutteme * prec;
703 return nutteme;
704
705 } // procedure truemean
706
707
708 // read data and apply corrections if needed from iau80rec data
709 // row 1-106, index 1-4
710 public static double rar80(int row, int index)
711 {
712 return iau80rec[row-1][index+4] * 0.0001 /3600.0;
713 }
714
715 // row 1-106, index 1-5
716 public static double iar80(int row, int index)
717 {
718 return iau80rec[row-1][index-1];
719 }
720
721 public static double[][] iau80rec =
722 {
723 {0,0,0,0,1,-171996,-174.2,92025,8.9,1}, // 1
724 {0,0,2,-2,2,-13187,-1.6,5736,-3.1,9},
725 {0,0,2,0,2,-2274,-0.2,977,-0.5,31},
726 {0,0,0,0,2,2062,0.2,-895,0.5,2},
727 {0,1,0,0,0,1426,-3.4,54,-0.1,10},
728 {1,0,0,0,0,712,0.1,-7,0,32},
729 {0,1,2,-2,2,-517,1.2,224,-0.6,11},
730 {0,0,2,0,1,-386,-0.4,200,0,33},
731 {1,0,2,0,2,-301,0,129,-0.1,34},
732 {0,-1,2,-2,2,217,-0.5,-95,0.3,12},
733 {1,0,0,-2,0,-158,0,-1,0,35},
734 {0,0,2,-2,1,129,0.1,-70,0,13},
735 {-1,0,2,0,2,123,0,-53,0,36},
736 {1,0,0,0,1,63,0.1,-33,0,38},
737 {0,0,0,2,0,63,0,-2,0,37},
738 {-1,0,2,2,2,-59,0,26,0,40},
739 {-1,0,0,0,1,-58,-0.1,32,0,39},
740 {1,0,2,0,1,-51,0,27,0,41},
741 {2,0,0,-2,0,48,0,1,0,14},
742 {-2,0,2,0,1,46,0,-24,0,3},
743 {0,0,2,2,2,-38,0,16,0,42},
744 {2,0,2,0,2,-31,0,13,0,45},
745 {2,0,0,0,0,29,0,-1,0,43},
746 {1,0,2,-2,2,29,0,-12,0,44},
747 {0,0,2,0,0,26,0,-1,0,46},
748 {0,0,2,-2,0,-22,0,0,0,15},
749 {-1,0,2,0,1,21,0,-10,0,47},
750 {0,2,0,0,0,17,-0.1,0,0,16},
751 {0,2,2,-2,2,-16,0.1,7,0,18},
752 {-1,0,0,2,1,16,0,-8,0,48},
753 {0,1,0,0,1,-15,0,9,0,17},
754 {1,0,0,-2,1,-13,0,7,0,49},
755 {0,-1,0,0,1,-12,0,6,0,19},
756 {2,0,-2,0,0,11,0,0,0,4},
757 {-1,0,2,2,1,-10,0,5,0,50},
758 {1,0,2,2,2,-8,0,3,0,54},
759 {0,-1,2,0,2,-7,0,3,0,53},
760 {0,0,2,2,1,-7,0,3,0,58},
761 {1,1,0,-2,0,-7,0,0,0,51},
762 {0,1,2,0,2,7,0,-3,0,52},
763 {-2,0,0,2,1,-6,0,3,0,20},
764 {0,0,0,2,1,-6,0,3,0,57},
765 {2,0,2,-2,2,6,0,-3,0,56},
766 {1,0,0,2,0,6,0,0,0,55},
767 {1,0,2,-2,1,6,0,-3,0,58},
768 {0,0,0,-2,1,-5,0,3,0,60},
769 {0,-1,2,-2,1,-5,0,3,0,21},
770 {2,0,2,0,1,-5,0,3,0,62},
771 {1,-1,0,0,0,5,0,0,0,61},
772 {1,0,0,-1,0,-4,0,0,0,24},
773 {0,0,0,1,0,-4,0,0,0,65},
774 {0,1,0,-2,0,-4,0,0,0,63},
775 {1,0,-2,0,0,4,0,0,0,64},
776 {2,0,0,-2,1,4,0,-2,0,22},
777 {0,1,2,-2,1,4,0,-2,0,23},
778 {1,1,0,0,0,-3,0,0,0,66},
779 {1,-1,0,-1,0,-3,0,0,0,6},
780 {-1,-1,2,2,2,-3,0,1,0,69},
781 {0,-1,2,2,2,-3,0,1,0,72},
782 {1,-1,2,0,2,-3,0,1,0,68},
783 {3,0,2,0,2,-3,0,1,0,71},
784 {-2,0,2,0,2,-3,0,1,0,5},
785 {1,0,2,0,0,3,0,0,0,67},
786 {-1,0,2,4,2,-2,0,1,0,82},
787 {1,0,0,0,2,-2,0,1,0,76},
788 {-1,0,2,-2,1,-2,0,1,0,74},
789 {0,-2,2,-2,1,-2,0,1,0,7},
790 {-2,0,0,0,1,-2,0,1,0,70},
791 {2,0,0,0,1,2,0,-1,0,75},
792 {3,0,0,0,0,2,0,0,0,77},
793 {1,1,2,0,2,2,0,-1,0,73},
794 {0,0,2,1,2,2,0,-1,0,78},
795 {1,0,0,2,1,-1,0,0,0,91},
796 {1,0,2,2,1,-1,0,1,0,85},
797 {1,1,0,-2,1,-1,0,0,0,102},
798 {0,1,0,2,0,-1,0,0,0,99},
799 {0,1,2,-2,0,-1,0,0,0,30},
800 {0,1,-2,2,0,-1,0,0,0,27},
801 {1,0,-2,2,0,-1,0,0,0,103},
802 {1,0,-2,-2,0,-1,0,0,0,100},
803 {1,0,2,-2,0,-1,0,0,0,94},
804 {1,0,0,-4,0,-1,0,0,0,80},
805 {2,0,0,-4,0,-1,0,0,0,83},
806 {0,0,2,4,2,-1,0,0,0,105},
807 {0,0,2,-1,2,-1,0,0,0,98},
808 {-2,0,2,4,2,-1,0,1,0,86},
809 {2,0,2,2,2,-1,0,0,0,90},
810 {0,-1,2,0,1,-1,0,0,0,101},
811 {0,0,-2,0,1,-1,0,0,0,97},
812 {0,0,4,-2,2,1,0,0,0,92},
813 {0,1,0,0,2,1,0,0,0,28},
814 {1,1,2,-2,2,1,0,-1,0,84},
815 {3,0,2,-2,2,1,0,0,0,93},
816 {-2,0,2,2,2,1,0,-1,0,81},
817 {-1,0,0,0,2,1,0,-1,0,79},
818 {0,0,-2,2,1,1,0,0,0,26},
819 {0,1,2,0,1,1,0,0,0,95},
820 {-1,0,4,0,2,1,0,0,0,87},
821 {2,1,0,-2,0,1,0,0,0,25},
822 {2,0,0,2,0,1,0,0,0,104},
823 {2,0,2,-2,1,1,0,-1,0,89},
824 {2,0,-2,0,1,1,0,0,0,8},
825 {1,-1,0,-2,0,1,0,0,0,88},
826 {-1,0,0,1,1,1,0,0,0,29},
827 {-1,-1,0,2,1,1,0,0,0,96},
828 {0,1,0,1,0,1,0,0,0,106}
829 };
830
831 }