Commit f9547e90 authored by Johan Richard's avatar Johan Richard
Browse files

Add reading RBF from .par

parent 1db83c27
......@@ -57,6 +57,7 @@ complex ci05f(double x, double y, double eps, double rc, double rcut);
complex ci10(double x, double y, double eps, double rc, double b0);
complex ci15(double x, double y, double eps, double rc, double b0);
void check_not_in_parallel(const char* s);
void checkpar();
double chiz(double z);
void classer(struct galaxie image[NFMAX][NIMAX], struct galaxie *cimage, long int *ni, long int *ncistart);
void cleanlens(double zimage);
......@@ -113,6 +114,7 @@ double d_random(int *idum);
double dratio(double zl, double zs);
void dratio_gal(struct galaxie *arclet, double zl);
double dratioprime(double zl, double zs);
void d_rbfsig(double **im, int nx, int ny, double scale);
double d_rndschechter(int *idum);
int d_rndtype(int *idum);
double d_rndz(double zmax, int *idum);
......
......@@ -107,6 +107,7 @@ struct g_observ
int setbin;
int bin;
int setseeing;
int setrbfsig;
int filtre;
double seeing;
double seeing_a;
......@@ -120,6 +121,7 @@ struct g_observ
double SKY;
double gain;
int idum;
double rbfsig;
};
struct g_mode
......
......@@ -9,13 +9,13 @@ AM_CFLAGS = @CFITSIO_HDR@ @WCS_HDR@ @GSL_HDR@ -I@top_srcdir@/include \
liblenstool_a_SOURCES = \
al_sq_point.c amplif.c amplif_mat.c amplif_matinv.c \
bayesapp.c bayesys3.c bicubic_coef.c bicubic_int.c \
brent.c chi_invim.c chsigne.c \
brent.c chi_invim.c chsigne.c \
classer.c cleanlens.c comp_chi_osv.c \
comp_dchi_osv.c copyright.c cor_seeing.c cor_shear.c \
cosmoratio.c cp_errdiff.c cp_images.c cp_im.c \
crea_filtre.c critic_an.c criticinv.c criticnew.c \
cv_cpsf.c d_binning.c d_bruiter.c d_gauss.c \
diag.c diff_mag.c d_integrer.c dist.c \
diag.c diff_mag.c d_integrer.c dist.c d_rbf.c \
dist_min.c distor.c do_itos.c d_poisson.c \
d_profil.c d_rndschechter.c d_rndtype.c d_rndz.c d_rndzsmail.c \
d_seeing.c e_amp.c ecrire_r.c e_dpl.c \
......
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include "lt.h"
/* Append seeing effect to the image im.
* Parameters :
* im : image in [ny][nx] format
* scale : pixel size in arcsec
*/
void d_rbf(double **im, int nx, int ny, double scale)
{
extern struct g_mode M;
const extern struct g_observ O;
register int i, j;
int n, N, err;
double **z, **tfr_z, **tfi_z, **tfr_f, **tfi_f; //f
double **tfr_conv, **tfi_conv; //**r_conv,**i_conv,
NPRINTF(stderr, "ADD: rbfseeing %lf\n",O.rbfsig);
/* on va inclure l'image dans une image plus grande de dimension n=2^N */
N = (int) ceil(log(((double)Max(nx, ny))) / log(2.));
n = (int) pow(2., ((double)N));
NPRINTF(stderr, "\tResize image in 2^N for FFT [%d,%d] --> [%d,%d]\n", nx, ny, n, n);
z = (double **) alloc_square_double(n, n); // z[ny][nx]
NPRINTF(stderr, "\tFill [%d,%d] image with original image\n", n, n);
for (i = 0; i < n; i++) // ny
for (j = 0; j < n; j++) // nx
{
if ((i < ny) && (j < nx))
z[i][j] = im[i][j];
else
z[i][j] = 0.;
}
NPRINTF(stderr, "\tAllocate images for Real & Imaginary part of the FT\n");
tfr_z = (double **) alloc_square_double(n, n);
tfi_z = (double **) alloc_square_double(n, n);
NPRINTF(stderr, "\tPerform the image FFT\n");
fftc_im(z, tfr_z, tfi_z, n, 1); // z is not used anymore for image
// By default O.filtre is null (see set_default.c) so do it!!!
if ( !O.filtre )
{
NPRINTF(stderr, "\tCreate the seeing filter\n");
printf("\tCreate the seeing filter\n");
if(O.setrbfsig==1) crea_filtre(O.rbfsig, scale, z, n); // z contains the seeing filter
}
//additional normalization of seeing filter
double psf_sum = 0;
for (i = 0 ; i < n ; i++)
for (j = 0 ; j < n ; j++)
psf_sum += z[i][j];
for (i = 0 ; i < n ; i++)
for (j = 0 ; j < n ; j++)
z[i][j] /= psf_sum;
tfr_f = (double **) alloc_square_double(n, n);
tfi_f = (double **) alloc_square_double(n, n);
NPRINTF(stderr, "\tPerform the filter FFT\n");
fftc_im(z, tfr_f, tfi_f, n, 1); // z is not used anymore for filter
free_square_double(z,n);
tfr_conv = (double **) alloc_square_double(n, n);
tfi_conv = (double **) alloc_square_double(n, n);
NPRINTF(stderr, "\tCompute the (image X filter) product in Fourier space\n");
ic_product(tfr_z, tfi_z, tfr_f, tfi_f, tfr_conv, tfi_conv, n, n);
// free_square_double(tfr_z,n);
// free_square_double(tfi_z,n);
free_square_double(tfr_f, n);
free_square_double(tfi_f, n);
//r_conv= (double **) alloc_square_double(n,n);
//i_conv= (double **) alloc_square_double(n,n);
NPRINTF(stderr, "\tCompute the IFFT of the product\n");
fftcc_im(tfr_conv, tfi_conv, tfr_z, tfi_z, n, -1);
free_square_double(tfr_conv, n);
free_square_double(tfi_conv, n);
NPRINTF(stderr, "\tSplit the image product of the convolution\n");
err = (int) split_image(tfr_z, n, n);
if (err != 0)
{
fprintf(stderr, "ERROR: when splitting the image\n");
exit(-1);
}
for (i = 0; i < ny; i++)
for (j = 0; j < nx; j++)
im[i][j] = tfr_z[i][j]; // [ny][nx]
free_square_double(tfr_z, n);
free_square_double(tfi_z, n);
// free_square_double(r_conv,n);
// free_square_double(i_conv,n);
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment