File: bias.cpp

package info (click to toggle)
pluto-find-orb 0.0~git20180227-2
  • links: PTS, VCS
  • area: main
  • in suites: experimental
  • size: 2,668 kB
  • sloc: cpp: 30,743; makefile: 263
file content (159 lines) | stat: -rw-r--r-- 6,193 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#include <stdlib.h>
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <assert.h>

void ra_dec_to_xy( const double ra, const double dec, double *x, double *y);
unsigned xy_to_healpix( const double x, const double y, const unsigned N);
int find_fcct_biases( const double ra, const double dec, const char catalog,
                 const double jd, double *bias_ra, double *bias_dec);

#define BIAS_NO_DATA_FOR_THAT_CATALOG     -1
#define BIAS_NO_BIAS_FILE                 -2
#define BIAS_NOT_ALLOCATED                -3
#define BIAS_FILE_WRONG_SIZE              -4

/* Compute star catalog position/proper motion biases,  from tables by
Farnocchia,  Chesley,  Chamberlin,  Tholen,  _Icarus_ 245 (2015) 94-111,
http://adsabs.harvard.edu/abs/2015Icar..245...94F  ,  "FCCT14".  The
actual debiasing file,  'bias.dat',  is in

ftp://ssd.jpl.nasa.gov/pub/ssd/debias/debias_2014.tgz

   The .tgz file contains two other files,  README.txt and tiles.dat.  As
far as Find_Orb is concerned,  you can ignore these.

   FCCT14 divides the sky into 12 * 64 * 64 = 49152 "tiles" of equal area and
passably square shape,  using the HEALPIX tesselation (see 'healpix.cpp').
Within each tile, the bias of 19 different catalogs in RA, dec,  and proper
motion in RA and dec are given.  (The biases are relative to a "presumed good
reference" subset of PPMXL.)  Those biases are stored in 'bias.dat'.

   This code reads in all 49152 * 19 biases (each of which is really
four biases,  one each in RA,  dec,  pm_RA,  pm_dec).  For a given
position,  it figures out the tile covering that point,  then extracts
the biases corresponding to the catalog in which we're interested and
computes the bias in RA and dec for a particular JD,  storing them in
*bias_ra and *bias_dec.       */

const char *fcct14_bias_file_name = NULL;
         /* override this if the bias file is elsewhere */

FILE *fopen_ext( const char *filename, const char *permits);   /* miscell.cpp */

int find_fcct_biases( const double ra, const double dec, const char catalog,
                 const double jd, double *bias_ra, double *bias_dec)
{
   static int16_t *bias_data = NULL;
   static const char *catalog_codes = "abcdegijlmopqruvwLN";
   const char *catalog_loc = strchr( catalog_codes, catalog);
         /* Number of HEALPIX tiles at level N=64 : */
   const int n_tiles = 12 * 64 * 64;
   const int n_cats = 19;
   int loc, i;
   double x, y, year;
   const double j2000 = 2451545.;      /* JD 2451545.0 = 2000 Jan 1.5 */
   static bool bias_file_known_to_be_missing = false;

   if( !bias_ra)         /* free up internal data */
      {
      if( bias_data)
         free( bias_data);
      bias_data = NULL;
      return( 0);
      }

   *bias_ra = *bias_dec = 0.;
   if( !catalog_loc)    /* don't have bias data for this catalog */
      return( BIAS_NO_DATA_FOR_THAT_CATALOG);
   if( bias_file_known_to_be_missing)
      return( BIAS_NO_BIAS_FILE);
   if( !bias_data)
      {
      FILE *ifile;
      const char *crunched_fcct14_file_name = "bias_bin.dat";
      size_t nbytes;
      const size_t csize = n_tiles * n_cats * 4 * sizeof( int16_t);
               /* The above "compressed size" corresponds to 12*64*64*19*4
                  = 3735552 short ints, or 7471104 bytes. */

      bias_data = (int16_t *)malloc( csize + 1);
      if( !bias_data)
         return( BIAS_NOT_ALLOCATED);
      ifile = fopen_ext( crunched_fcct14_file_name, "crb");
      if( ifile)     /* not our first time:  we have a compressed file */
         {
         nbytes = fread( bias_data, 1, csize + 1, ifile);
         fclose( ifile);
         if( nbytes != csize)
            return( BIAS_FILE_WRONG_SIZE);
         }
      else
         {
         char buff[600];   /* first time: read ASCII file & binary-ize; */
         FILE *ofile;     /* store binary version for all future use   */

         ifile = NULL;
         if( fcct14_bias_file_name && * fcct14_bias_file_name)
            ifile = fopen( fcct14_bias_file_name, "rb");
         if( !ifile)
            ifile = fopen_ext( "bias.dat", "crb");
         if( !ifile)
            {
            bias_file_known_to_be_missing = true;
            free( bias_data);
            return( BIAS_NO_BIAS_FILE);
            }
         loc = 0;
         while( loc < n_tiles && fgets( buff, sizeof( buff), ifile))
            if( *buff != '!')
               {
               char *tptr = buff;
               int16_t *bptr = bias_data + loc * n_cats * 4;

               for( i = 0; i < n_cats * 4; i++, bptr++)
                  {
                  bool is_negative = false;

                  while( *tptr == ' ')
                     tptr++;
                  if( *tptr == '-')
                     {
                     is_negative = true;
                     tptr++;
                     }
                  *bptr = 0;
                  while( *tptr != ' ')
                     {
                     if( *tptr >= '0' && *tptr <= '9')
                        *bptr = (int16_t)( *bptr * 10 + *tptr - '0');
                     tptr++;
                     }
                  if( is_negative)
                     *bptr = -*bptr;
                  }
               loc++;
               }
         fclose( ifile);
         assert( loc == n_tiles);
         ofile = fopen_ext( crunched_fcct14_file_name, "fcwb");
         nbytes = fwrite( bias_data, 1, csize, ofile);
         assert( nbytes == csize);
         fclose( ofile);
         }
      }
   ra_dec_to_xy( ra, dec, &x, &y);
   loc = xy_to_healpix( x, y, 64);
// debug_printf( "catalog %c, tile %d: ", catalog, loc);
   loc = loc * n_cats * 4 + (int)(catalog_loc - catalog_codes) *4;
   year = (jd - j2000) / 365.25;
// debug_printf( "%d %d %d %d\n", bias_data[loc], bias_data[loc + 1],
//                   bias_data[loc + 2], bias_data[loc + 3]);
   *bias_ra  =  (double)bias_data[loc]     / 1000.     /* RA bias, in arcsec */
       + year * (double)bias_data[loc + 2] / 100000.;  /* RA pm bias, in mas/yr */
   *bias_dec =  (double)bias_data[loc + 1] / 1000.     /* dec bias, in arcsec */
       + year * (double)bias_data[loc + 3] / 100000.;  /* dec pm bias, in mas/yr */
   return( 0);
}