o_pixsource.c 7.42 KB
Newer Older
1
#include<stdio.h>
2
#include<stdlib.h>
3
4
5
6
7
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
8
#include<lt.h>
9
10
11
12
13
14
15
/****************************************************************/
/*      nom:        o_pixsource              */
/*      auteur:     Johan Richard         */
/*      date:       02/01/14            */
/*      place:      Lyon                */
/****************************************************************/

16
void precompsource(){
17
18
    FILE *fp;
    struct point Pi, Ps;
19
   const extern  struct  pot lens[];
20
    const extern struct g_source  S;
21
22
23
24
   extern struct g_mode M; // M.iclean
   const extern struct g_observ  O;
   extern struct g_pixel imFrame,ps,PSF;
   extern double **fluxmap;
Johan Richard's avatar
Johan Richard committed
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
    double **tmpimage,xx ;
  extern   int **nimage;
  const extern struct g_frame   F;  

  // Setting the imframe parameters    
/*  imFrame.xmin = F.xmin;
    imFrame.xmax = F.xmax;
    imFrame.ymin = F.ymin;
    imFrame.ymax = F.ymax;
    imFrame.nx = M.npixel;   
    double    scale = (imFrame.xmax-imFrame.xmin) / (imFrame.nx-1.);
    imFrame.ny = round((imFrame.ymax-imFrame.ymin )/ scale + 1.);  // to prevent trunc(ny) behavior
    imFrame.pixelx = scale;
    imFrame.pixely = scale;
    imFrame.ech = 10;  */
    xx = (double)imFrame.ech; 

//Setting the ps parameters
double dlsds = dratio(lens[0].z, M.zclean);
   fprintf(stderr,"dlsds=%lf\n",dlsds);
 if( strcmp(M.centerfile, "") )
    {
        ps.pixelx = imFrame.pixelx / ps.ech;
        ps.pixely = imFrame.pixelx/ ps.ech;
        s_sourcebox(&ps, M.centerfile, dlsds, M.zclean);
    } 

52
53
   int ni=imFrame.nx*imFrame.ny;
   int ns=ps.nx*ps.ny;
54
    double stepx,stepxs,stepy,stepys;
55
   fluxmap=(double **)alloc_square_double(ni,ns);
Johan Richard's avatar
Johan Richard committed
56
   
57
   double sizeratio=imFrame.pixelx/ps.pixelx; // pixel size ratio between image and source plane
58
  double xs,ys;
Johan Richard's avatar
Johan Richard committed
59
60
61
62
     fprintf(stderr,"ps.xmax=%lf\t ps.ymax=%lf\n",ps.xmax,ps.ymax);
    fprintf(stderr,"ps.xmin=%lf\t ps.ymin=%lf\n",ps.xmin,ps.ymin);
   fprintf(stderr,"imframe.nx=%d\t imFrame.ny=%d\n",imFrame.nx,imFrame.ny);
   fprintf(stderr,"imFrame.pixelx=%lf\t imFrame.pixely=%lf\n",imFrame.pixelx,imFrame.pixely);
63
64
65
66
    fprintf(stderr,"ps.pixelx=%lf\t ps.pixely=%lf\n",ps.pixelx,ps.pixely);
    fprintf(stderr,"imFrame.xmax=%lf\t imFrame.ymax=%lf\n",imFrame.xmax,imFrame.ymax);
    fprintf(stderr,"imFrame.xmin=%lf\t imFrame.ymin=%lf\n",imFrame.xmin,imFrame.ymin);
    fprintf(stderr, "sizeratio= %lf\n", sizeratio);
Johan Richard's avatar
Johan Richard committed
67
    fprintf(stderr,"ech=%d\n",imFrame.ech);
68
   tmpimage=(double **)alloc_square_double(imFrame.ny,imFrame.nx);
Johan Richard's avatar
Johan Richard committed
69
70
    nimage=(int **)alloc_square_int(imFrame.ny,imFrame.nx);
   int is,js,i,j,ii,jj;
71
72
73
74
75
76
77
78
79
    double z;
    stepx=(imFrame.xmax-imFrame.xmin)/(imFrame.nx-1);
    stepy=(imFrame.ymax-imFrame.ymin)/(imFrame.ny-1);
    stepxs=(ps.xmax-ps.xmin)/(ps.nx-1);
    stepys=(ps.ymax-ps.ymin)/(ps.ny-1);
    fprintf(stderr,"stepxs=%.8lf\t stepys=%.8lf\n",stepxs,stepys);
    fprintf(stderr,"stepx=%.8lf\t stepy=%.8lf\n",stepx,stepy);
    for (i=0;i<imFrame.nx;i++)
    for(j=0;j<imFrame.ny;j++)
Johan Richard's avatar
Johan Richard committed
80
       {
81
        tmpimage[j][i]=0.0;
Johan Richard's avatar
Johan Richard committed
82
83
84
       nimage[j][i]=0;
       }
        
85
86
87
88
89
    for (i=0;i<ni;i++)
        for(j=0;j<ns;j++)
            fluxmap[i][j]=0.0;
    xs= ps.xmin;
  //  ys=ps.ymin;
Johan Richard's avatar
Johan Richard committed
90
//    fp = fopen("fluxmap.dat","w");
91
92
93
     for(is=0;is<ps.nx;is++)
   { fprintf(stderr, "is=%d\n",is);
       ys=ps.ymin;
94
      for(js=0;js<ps.ny;js++)
Johan Richard's avatar
Johan Richard committed
95
96
97
98
99
100
101
      {  //Pi.x= imFrame.xmin - stepx/2;
       Pi.x= imFrame.xmin - stepx/2 * (1.- 1./xx);
          for (i=0;i<imFrame.nx;i++)
           for(j=0;j<imFrame.ny;j++)
       {
        tmpimage[j][i]=0.0;
       }   
102
103
        for(i=0;i<imFrame.nx;i++)
        {
Johan Richard's avatar
Johan Richard committed
104
105
106
         for (ii = 0; ii < imFrame.ech; ii++)
         { //  Pi.y=imFrame.ymin - stepy/2;
        Pi.y=imFrame.ymin - stepy/2 * (1.- 1./xx); 
107
108
          for(j=0;j<imFrame.ny;j++)
          {
Johan Richard's avatar
Johan Richard committed
109
110
                for (jj = 0; jj < imFrame.ech; jj++)
                {
111
            e_dpl(&Pi, dlsds,S.zs, &Ps);
Johan Richard's avatar
Johan Richard committed
112
113
114
115
116
117
118
119
120
                 z= fmax(fabs(Ps.x-xs)/stepxs,fabs(Ps.y-ys)/stepys);
              if(z<0.5)
              { 
                   nimage[j][i]++;  
               //   fprintf(stderr,"nimage[%d][%d]=%d\n",j,i,nimage[j][i]);
                  tmpimage[j][i]+=1.0;//(sizeratio*sizeratio);
               } 
              Pi.y +=stepy/xx;
                 }
121
          }
Johan Richard's avatar
Johan Richard committed
122
123
124
125
            Pi.x +=stepx/xx;
         }
        }  
      /*  convolve tmpimage                                  with PSF */
126
127
        if (O.setseeing)
               d_seeing_omp(tmpimage, imFrame.nx, imFrame.ny, imFrame.pixelx);
128
       /* store tmpimage in flux matrix */
129
        for(i=0;i<imFrame.nx;i++)
Johan Richard's avatar
Johan Richard committed
130
        { 
131
          for(j=0;j<imFrame.ny;j++)
Johan Richard's avatar
Johan Richard committed
132
133
134
135
          {  
          fluxmap[i*imFrame.ny+j][is*ps.ny+js]=tmpimage[j][i]/xx/xx;
   //   if(fluxmap[i*imFrame.ny+j][is*ps.ny+js]>1.e-15)
      //     fprintf(fp,"fluxmap[%d][%d]=%e\t tmpimage[%d][%d]=%e\n",(i*imFrame.ny+j),(is*ps.ny+js),fluxmap[i*imFrame.ny+j][is*ps.ny+js],j,i,tmpimage[j][i]);
136
          }
137
        }
138
        ys+=stepys;
139
     }
140
       xs+=stepxs;
141
   }
Johan Richard's avatar
Johan Richard committed
142
143
144
145
  //    fclose(fp);

  wri_fits_abs("npixel_abs.fits", nimage, imFrame.nx, imFrame.ny, imFrame.xmin, imFrame.xmax, imFrame.ymin, imFrame.ymax, M.ref_ra, M.ref_dec);
  
146
   free_square_double(tmpimage,imFrame.ny);
147
148
149
   
}

Johan Richard's avatar
Johan Richard committed
150
void testpixel(char *imname,char *sname,double *sflux)
151
152
153
154
155
156
{
    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;
Johan Richard's avatar
Johan Richard committed
157
158
    double **zz,**source;
    double f, scale,xx;
159
160
161
    int nx, ny;
    FILE *fp;
    int i,j,k;
Johan Richard's avatar
Johan Richard committed
162
     extern   int **nimage;
163
164
165
166
167
168
169
double  ra, dec, width, height;
    xmin = imFrame.xmin;
    xmax = imFrame.xmax;
    ymin = imFrame.ymin;
    ymax = imFrame.ymax;

    dx = xmax - xmin;
Johan Richard's avatar
Johan Richard committed
170
   dy = ymax - ymin;
171
    nx = imFrame.nx;// default: assign np to nx and adjust ny
Johan Richard's avatar
Johan Richard committed
172
173
174
  //  ny = imFrame.ny;
   scale = dx / (nx-1);
 ny =  round(dy / scale + 1.);  // to prevent trunc(ny) behaviour
175
176
177
178
179
180
181
182
183
184

    dx = (nx - 1) * scale;
    dy = (ny - 1) * scale;
    xmax = xmin + dx;
    ymax = ymin + dy;
    int ni=imFrame.nx*imFrame.ny;
    int ns=ps.nx*ps.ny;
  fprintf(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);
Johan Richard's avatar
Johan Richard committed
185
186
187
188
189
  fprintf(stderr,"nx=%d\t ny=%d\n",nx,ny);
  char filename[128];
  sprintf(filename, "params_%s.dat", imname);
 //   fprintf(stderr,"filename = %s\n",filename);
  //  fp = fopen(filename,"w");
190
    for(i=0;i<imFrame.nx;i++)
Johan Richard's avatar
Johan Richard committed
191
    { 
192
193
       for(j=0;j<imFrame.ny;j++)
       {  zz[j][i]=0.0;
Johan Richard's avatar
Johan Richard committed
194
     //  xx = (double)nimage[j][i];
195
          for(k=0;k<ns;k++)
Johan Richard's avatar
Johan Richard committed
196
         { //fprintf(fp,"fluxmap[%d][%d]=%e\t sflux[%d]=%e\t",i*ny+j,k,fluxmap[i*ny+j][k],k,sflux[k]);
197
          zz[j][i]+=fluxmap[i*ny+j][k]*sflux[k];
Johan Richard's avatar
Johan Richard committed
198
199
200
201
202
          }
        //   fprintf(fp,"zz[%d][%d]=%e\t fluxmap[%d][%d]=%e\n",j,i,zz[j][i]);
      //    if (nimage[j][i]>1)
        // zz[j][i]=zz[j][i]/xx;
          
203
204
       }
    }
Johan Richard's avatar
Johan Richard committed
205
206
207
    
  //  fclose(fp);
     wrf_fits_abs(imname, zz, nx, ny, xmin, xmax, ymin, ymax, M.ref_ra, M.ref_dec);
208
   
Johan Richard's avatar
Johan Richard committed
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
    source = (double **) alloc_square_double(ps.ny, ps.nx);
    // Also write the corresponding source file using the ps parameters
 
    for(i=0;i<ps.nx;i++)
    {
        for(j=0;j<ps.ny;j++)
        
  source[j][i]=sflux[i*ps.ny+j];
        
    }
 //   Nprintf(stderr,"source[%d][%d]=%e\n",9,19,source[9][19]);
   wrf_fits_abs(sname, source, ps.nx, ps.ny, ps.xmin, ps.xmax, ps.ymin, ps.ymax, M.ref_ra, M.ref_dec);
    
    free_square_double(zz, ny);
    free_square_double(source, ps.ny);
224
}