o_pixsource.c 3.41 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
#include<stdio.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>


/****************************************************************/
/*      nom:        o_pixsource              */
/*      auteur:     Johan Richard         */
/*      date:       02/01/14            */
/*      place:      Lyon                */
/****************************************************************/

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
void precompsource(){
   const extern  struct  pot lens[];
   extern struct g_mode M; // M.iclean
   const extern struct g_observ  O;
   extern struct g_pixel imFrame,ps,PSF;
   extern double **fluxmap;

   int ni=imFrame.nx*imFrame.ny;
   int ns=ps.nx*ps.ny;

   fluxmap=(double **)alloc_square_double(ni,ns);

   double dlsds = dratio(lens[0].z, M.zclean);
   double sizeratio=imFrame.pixelx/ps.pixelx; // pixel size ratio between image and source plane
   double xi,yi,xs,ys,xs2,ys2;

   double **tmpimage=(double **)alloc_square_double(imFrame.nx,imFrame.ny);
   int is,js,i,j;

   for(is=0;is<ps.nx;is++)
   {
      for(js=0;js<ps.ny;js++)
      {
        xs = is * (ps.xmax - ps.xmin) / (ps.nx - 1) + ps.xmin;
        ys = js * (ps.ymax - ps.ymin) / (ps.ny - 1) + ps.ymin;
        for(i=0;i<imFrame.nx;i++)
        {
          for(j=0;j<imFrame.ny;j++)
          {
            xi = i * (imFrame.xmax - imFrame.xmin) / (imFrame.nx - 1) + imFrame.xmin;
            yi = j * (imFrame.ymax - imFrame.ymin) / (imFrame.ny - 1) + imFrame.ymin;
            //e_dpl(xi,yi, dlsds, xs2,ys2);
            if((xs2-xs)*(xs2-xs)+(ys2-ys)*(ys2-ys)<0.25) tmpimage[i][j]=sizeratio;
            else tmpimage[i][j]=0.0;
          }
        }
        //convolve tmpimage with PSF
        if (O.setseeing)
               d_seeing_omp(tmpimage, imFrame.nx, imFrame.ny, imFrame.pixelx);

        //store tmpimage in flux matrix
        for(i=0;i<imFrame.nx;i++)
        {
          for(j=0;j<imFrame.ny;j++)
          {
              fluxmap[i*imFrame.ny+j][is*ps.ny+js]=tmpimage[i][j];
          }   
        }
     }
   }
   free_square_double(tmpimage,imFrame.ny);
}

void    testpixel(double *sflux)
{
    const extern struct g_observ  O;
    extern struct g_mode    M;
    extern struct g_pixel   ps,imFrame;
    extern double **fluxmap;

    double xmin, ymin, xmax, ymax, dx, dy;
    double f, scale;
    int nx, ny, ny_sav;
    double **zz;
    int i,j,k;

    xmin = imFrame.xmin;
    xmax = imFrame.xmax;
    ymin = imFrame.ymin;
    ymax = imFrame.ymax;

    dx = xmax - xmin;
    dy = ymax - ymin;
    nx = imFrame.nx;    // default: assign np to nx and adjust ny
    scale = dx / (nx-1);
    ny = dy / scale + 1;  // warning: trunc(ny) behavior 

    dx = (nx - 1) * scale;
    dy = (ny - 1) * scale;
    xmax = xmin + dx;
    ymax = ymin + dy;

    NPRINTF(stderr, "\timage (%d,%d) s=%.3lf [%.3lf,%.3lf] [%.3lf,%.3lf]\n",
            nx, ny, scale, xmin, xmax, ymin, ymax);

    zz = (double **) alloc_square_double(ny, nx);

    int ns=ps.nx*ps.ny;
    for(i=0;i<nx;i++)
    {
       for(j=0;j<ny;j++)
       {
          zz[i][j]=0.0;
          for(k=0;k<ns;k++)
          {
             zz[i][j]+=fluxmap[k][i*ny+j]*sflux[k];
          }
       }
    }

    if (O.setseeing)
        d_seeing(zz, nx, ny, scale);

    if ( M.iref > 0 )
        wrf_fits_abs("testpixel.fits", zz, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
    else
        wrf_fits("testpixel.fits", zz, nx, ny, xmin, xmax, ymin, ymax);

    free_square_double(zz, ny_sav);

126
}