IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37674


Ignore:
Timestamp:
Nov 26, 2014, 6:21:42 AM (12 years ago)
Author:
eugene
Message:

use double for gaussj inversion

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20140904/Ohana/src/libdvo/src/AstromOffsetMapOps.c

    r37445 r37674  
    33/* AstromOffsetMap functions:
    44 */
     5
     6int dump_map_data (float *x, float *y, float *f, int Npts, char *filename);
    57
    68float AstromOffsetMapValue (AstromOffsetMap *map, float x, float y, int xdir) {
    … …  
    108110  int Ny = map->Ny;
    109111
    110   float **A, **B;
    111   ALLOCATE (A, float *, Nx*Ny);
    112   ALLOCATE (B, float *, Nx*Ny);
     112  double **A, **B;
     113  ALLOCATE (A, double *, Nx*Ny);
     114  ALLOCATE (B, double *, Nx*Ny);
    113115  for (i = 0; i < Nx*Ny; i++) {
    114     ALLOCATE (A[i], float, Nx*Ny);
    115     memset (A[i], 0, sizeof(float)*Nx*Ny);
    116     ALLOCATE (B[i], float, 1);
    117     memset (B[i], 0, sizeof(float));
     116    ALLOCATE (A[i], double, Nx*Ny);
     117    memset (A[i], 0, sizeof(double)*Nx*Ny);
     118    ALLOCATE (B[i], double, 1);
     119    memset (B[i], 0, sizeof(double));
    118120  }   
    119121
    … …  
    135137    for (iy = 0; iy < Ny; iy++) {
    136138      // define & init summing variables
    137       float rx_rx_ry_ry = 0;
    138       float rx_rx_dy_ry = 0;
    139       float dx_rx_ry_ry = 0;
    140       float dx_rx_dy_ry = 0;
    141       float fi_rx_ry    = 0;
    142       float rx_rx_py_py = 0;
    143       float rx_rx_qy_py = 0;
    144       float dx_rx_py_py = 0;
    145       float dx_rx_qy_py = 0;
    146       float fi_rx_py    = 0;
    147       float px_px_ry_ry = 0;
    148       float px_px_dy_ry = 0;
    149       float qx_px_ry_ry = 0;
    150       float qx_px_dy_ry = 0;
    151       float fi_px_ry    = 0;
    152       float px_px_py_py = 0;
    153       float px_px_qy_py = 0;
    154       float qx_px_py_py = 0;
    155       float qx_px_qy_py = 0;
    156       float fi_px_py    = 0;
     139      double rx_rx_ry_ry = 0;
     140      double rx_rx_dy_ry = 0;
     141      double dx_rx_ry_ry = 0;
     142      double dx_rx_dy_ry = 0;
     143      double fi_rx_ry    = 0;
     144      double rx_rx_py_py = 0;
     145      double rx_rx_qy_py = 0;
     146      double dx_rx_py_py = 0;
     147      double dx_rx_qy_py = 0;
     148      double fi_rx_py    = 0;
     149      double px_px_ry_ry = 0;
     150      double px_px_dy_ry = 0;
     151      double qx_px_ry_ry = 0;
     152      double qx_px_dy_ry = 0;
     153      double fi_px_ry    = 0;
     154      double px_px_py_py = 0;
     155      double px_px_qy_py = 0;
     156      double qx_px_py_py = 0;
     157      double qx_px_qy_py = 0;
     158      double fi_px_py    = 0;
    157159
    158160      // generate the sums for the fitting matrix element I,J
    … …  
    167169
    168170        // base coordinate offset for this point (x,y) relative to this map element (n,m)
    169         // float dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
    170         // float dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
    171 
    172         float dx = x[i] * map->dX - ix - 0.5;
    173         float dy = y[i] * map->dY - iy - 0.5;
     171        // double dx = psImageBinningGetRuffX (map->binning, x->data.F32[i]) - (n + 0.5);
     172        // double dy = psImageBinningGetRuffY (map->binning, y->data.F32[i]) - (m + 0.5);
     173
     174        double dx = x[i] * map->dX - ix - 0.5;
     175        double dy = y[i] * map->dY - iy - 0.5;
    174176
    175177        // XXX do I need to do something different at the edge or not?  I want to allow
    … …  
    195197
    196198        // related offset values
    197         float rx = 1.0 - dx;
    198         float ry = 1.0 - dy;
    199         float px = 1.0 + dx;
    200         float py = 1.0 + dy;
    201         float qx = 0.0 - dx;
    202         float qy = 0.0 - dy;
     199        double rx = 1.0 - dx;
     200        double ry = 1.0 - dy;
     201        double px = 1.0 + dx;
     202        double py = 1.0 + dy;
     203        double qx = 0.0 - dx;
     204        double qy = 0.0 - dy;
    203205
    204206        // sum the appropriate elements for the different quadrants
    … …  
    306308  }
    307309
    308   if (0) {
     310  if (1) {
    309311    FILE *fd = fopen ("matrix.dat", "w");
    310312    for (i = 0; i < Nx*Ny; i++) {
    … …  
    317319  }
    318320
    319   if (!fgaussjordan(A, B, Nx*Ny, 1)) {
     321  if (!dgaussjordan(A, B, Nx*Ny, 1)) {
    320322    fprintf (stderr, "FAIL\n");
    321323    exit (1);
    … …  
    350352}
    351353
     354int dump_map_data (float *x, float *y, float *f, int Npts, char *filename) {
     355
     356  FILE *fout = fopen (filename, "w");
     357 
     358  int i;
     359  for (i = 0; i < Npts; i++) {
     360    fprintf (fout, "%d %f %f %f\n", i, x[i], y[i], f[i]);
     361  }
     362  fclose (fout);
     363  return TRUE;
     364}
     365
Note: See TracChangeset for help on using the changeset viewer.