Commit aa3c2f14 authored by Eric Jullo's avatar Eric Jullo
Browse files

Add wcs support for galaxy catalog and healpix shear maps

parent e9db396c
# generated automatically by aclocal 1.15 -*- Autoconf -*-
# generated automatically by aclocal 1.15.1 -*- Autoconf -*-
# Copyright (C) 1996-2014 Free Software Foundation, Inc.
# Copyright (C) 1996-2017 Free Software Foundation, Inc.
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -20,7 +20,7 @@ You have another version of autoconf. It may work, but is not guaranteed to.
If you have problems, you may need to regenerate the build system entirely.
To do so, use the procedure documented by the package, typically 'autoreconf'.])])
# Copyright (C) 2002-2014 Free Software Foundation, Inc.
# Copyright (C) 2002-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -35,7 +35,7 @@ AC_DEFUN([AM_AUTOMAKE_VERSION],
[am__api_version='1.15'
dnl Some users find AM_AUTOMAKE_VERSION and mistake it for a way to
dnl require some minimum version. Point them to the right macro.
m4_if([$1], [1.15], [],
m4_if([$1], [1.15.1], [],
[AC_FATAL([Do not call $0, use AM_INIT_AUTOMAKE([$1]).])])dnl
])
......@@ -51,14 +51,14 @@ m4_define([_AM_AUTOCONF_VERSION], [])
# Call AM_AUTOMAKE_VERSION and AM_AUTOMAKE_VERSION so they can be traced.
# This function is AC_REQUIREd by AM_INIT_AUTOMAKE.
AC_DEFUN([AM_SET_CURRENT_AUTOMAKE_VERSION],
[AM_AUTOMAKE_VERSION([1.15])dnl
[AM_AUTOMAKE_VERSION([1.15.1])dnl
m4_ifndef([AC_AUTOCONF_VERSION],
[m4_copy([m4_PACKAGE_VERSION], [AC_AUTOCONF_VERSION])])dnl
_AM_AUTOCONF_VERSION(m4_defn([AC_AUTOCONF_VERSION]))])
# AM_AUX_DIR_EXPAND -*- Autoconf -*-
# Copyright (C) 2001-2014 Free Software Foundation, Inc.
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -110,7 +110,7 @@ am_aux_dir=`cd "$ac_aux_dir" && pwd`
# AM_CONDITIONAL -*- Autoconf -*-
# Copyright (C) 1997-2014 Free Software Foundation, Inc.
# Copyright (C) 1997-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -141,7 +141,7 @@ AC_CONFIG_COMMANDS_PRE(
Usually this means the macro was only invoked conditionally.]])
fi])])
# Copyright (C) 1999-2014 Free Software Foundation, Inc.
# Copyright (C) 1999-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -332,7 +332,7 @@ _AM_SUBST_NOTMAKE([am__nodep])dnl
# Generate code to set up dependency tracking. -*- Autoconf -*-
# Copyright (C) 1999-2014 Free Software Foundation, Inc.
# Copyright (C) 1999-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -408,7 +408,7 @@ AC_DEFUN([AM_OUTPUT_DEPENDENCY_COMMANDS],
# Do all the work for Automake. -*- Autoconf -*-
# Copyright (C) 1996-2014 Free Software Foundation, Inc.
# Copyright (C) 1996-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -605,7 +605,7 @@ for _am_header in $config_headers :; do
done
echo "timestamp for $_am_arg" >`AS_DIRNAME(["$_am_arg"])`/stamp-h[]$_am_stamp_count])
# Copyright (C) 2001-2014 Free Software Foundation, Inc.
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -626,7 +626,7 @@ if test x"${install_sh+set}" != xset; then
fi
AC_SUBST([install_sh])])
# Copyright (C) 2003-2014 Free Software Foundation, Inc.
# Copyright (C) 2003-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -647,7 +647,7 @@ AC_SUBST([am__leading_dot])])
# Check to see how 'make' treats includes. -*- Autoconf -*-
# Copyright (C) 2001-2014 Free Software Foundation, Inc.
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -697,7 +697,7 @@ rm -f confinc confmf
# Fake the existence of programs that GNU maintainers use. -*- Autoconf -*-
# Copyright (C) 1997-2014 Free Software Foundation, Inc.
# Copyright (C) 1997-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -736,7 +736,7 @@ fi
# Helper functions for option handling. -*- Autoconf -*-
# Copyright (C) 2001-2014 Free Software Foundation, Inc.
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -765,7 +765,7 @@ AC_DEFUN([_AM_SET_OPTIONS],
AC_DEFUN([_AM_IF_OPTION],
[m4_ifset(_AM_MANGLE_OPTION([$1]), [$2], [$3])])
# Copyright (C) 1999-2014 Free Software Foundation, Inc.
# Copyright (C) 1999-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -812,7 +812,7 @@ AC_LANG_POP([C])])
# For backward compatibility.
AC_DEFUN_ONCE([AM_PROG_CC_C_O], [AC_REQUIRE([AC_PROG_CC])])
# Copyright (C) 2001-2014 Free Software Foundation, Inc.
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -831,7 +831,7 @@ AC_DEFUN([AM_RUN_LOG],
# Check to make sure that the build environment is sane. -*- Autoconf -*-
# Copyright (C) 1996-2014 Free Software Foundation, Inc.
# Copyright (C) 1996-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -912,7 +912,7 @@ AC_CONFIG_COMMANDS_PRE(
rm -f conftest.file
])
# Copyright (C) 2009-2014 Free Software Foundation, Inc.
# Copyright (C) 2009-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -972,7 +972,7 @@ AC_SUBST([AM_BACKSLASH])dnl
_AM_SUBST_NOTMAKE([AM_BACKSLASH])dnl
])
# Copyright (C) 2001-2014 Free Software Foundation, Inc.
# Copyright (C) 2001-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -1000,7 +1000,7 @@ fi
INSTALL_STRIP_PROGRAM="\$(install_sh) -c -s"
AC_SUBST([INSTALL_STRIP_PROGRAM])])
# Copyright (C) 2006-2014 Free Software Foundation, Inc.
# Copyright (C) 2006-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......@@ -1019,7 +1019,7 @@ AC_DEFUN([AM_SUBST_NOTMAKE], [_AM_SUBST_NOTMAKE($@)])
# Check how to create a tarball. -*- Autoconf -*-
# Copyright (C) 2004-2014 Free Software Foundation, Inc.
# Copyright (C) 2004-2017 Free Software Foundation, Inc.
#
# This file is free software; the Free Software Foundation
# gives unlimited permission to copy and/or distribute it,
......
......@@ -14,7 +14,7 @@
#define NMAX 5000 // maximum number of segments for the critical lines
#define NPMAX 5000
#define NPZMAX 9 /* maximum number of critical lines in g_cline struct*/
#define NLMAX 4000 // maximum number of clumps in the lens[] array
#define NLMAX 70000 // maximum number of clumps in the lens[] array
#define NIMAX 50 /* maximum images per family */
#define NFMAX 908 /* maximum number of families */
#define NPAMAX 35 // number of free parameters (see #define in structure.h)
......
#include "wcs.h"
/*********************************************************************
* Simple interpolator of a wcs transformation
*
* Author: Eric Jullo
* Date: July 2019
*
********************************************************************/
#define NRR 100 // Number of sampling points to initialize the splines
/********************************************************************
* Initialize a GSL interp spline interpolation
*/
int wcs_alloc(struct WorldCoor *wcsinfo, double ramin, double ramax,
double decmin, double decmax, int nx, int ny)
{
gsl_interp_accel *acc1; // ra interp
gsl_spline *hf1;
gsl_interp_accel *acc2; // dec interp
gsl_spline *hf2;
int i;
double arr1[NRR], arr2[NRR], x, y;
const gsl_interp2d_type *T = gsl_interp2d_bilinear;
double *ra = malloc(nx * ny * sizeof(double));
gsl_spline2d *spline = gsl_spline2d_alloc(T, nx, ny);
acc1=gsl_interp2d_accel_alloc();
hf1=gsl_spline_alloc(gsl_interp_cspline,nx);
acc2=gsl_interp_accel_alloc();
hf2=gsl_spline_alloc(gsl_interp_cspline,ny);
for (i=0; i<NRR; i++) {
pix2wcs(wcsinfo, i, j
}
}
......@@ -75,23 +75,23 @@ static void cpa2d_free(struct cpa2d* c)
free(c);
}
//
inline void cpa2d_set(struct cpa2d* c, int i, int j, double re, double im)
void cpa2d_set(struct cpa2d* c, int i, int j, double re, double im)
{
c->v[(i * c->ny + j) * 2] = re;
c->v[(i * c->ny + j) * 2 + 1] = im;
}
//
inline void cpa2d_set_re(struct cpa2d* c, int i, int j, double re)
void cpa2d_set_re(struct cpa2d* c, int i, int j, double re)
{
c->v[(i * c->ny + j) * 2] = re;
}
//
inline double cpa2d_get_re(struct cpa2d* c, int i, int j)
double cpa2d_get_re(struct cpa2d* c, int i, int j)
{
return c->v[(i * c->ny + j) * 2];
}
//
inline double cpa2d_get_im(struct cpa2d* c, int i, int j)
double cpa2d_get_im(struct cpa2d* c, int i, int j)
{
return c->v[(i * c->ny + j) * 2 + 1];
}
......
......@@ -850,15 +850,15 @@ if ( ilens->epot > 2E-25 )
g2.c = kap + gamma.x;
break;
case(122): // NFW + large-scale structure (correlated linear spectrum)
Qe.x = sqrt(1. - ilens->emass) * Q.x;
Qe.y = sqrt(1. + ilens->emass) * Q.y;
QP = polxy(Qe);
if ( QP.r == 0 ) QP.r = 1E-8;
kap = nfw_kappa_eps(QP.r, ilens->rc, QP.theta, ilens->b0, ilens->emass);
gamma = nfw_gamma_eps(QP.r, ilens->rc, QP.theta, ilens->b0, ilens->emass);
g2.a = kap - gamma.x;
g2.b = g2.d = gamma.y;
g2.c = kap + gamma.x;
//Qe.x = sqrt(1. - ilens->emass) * Q.x;
//Qe.y = sqrt(1. + ilens->emass) * Q.y;
//QP = polxy(Qe);
//if ( QP.r == 0 ) QP.r = 1E-8;
//kap = nfw_kappa_eps(QP.r, ilens->rc, QP.theta, ilens->b0, ilens->emass);
//gamma = nfw_gamma_eps(QP.r, ilens->rc, QP.theta, ilens->b0, ilens->emass);
//g2.a = kap - gamma.x;
//g2.b = g2.d = gamma.y;
//g2.c = kap + gamma.x;
QP = polxy(Q);
const extern struct g_cosmo C;
......
......@@ -167,7 +167,7 @@ void f_shape( long int *istart,
liste[i].C.x -= M.ref_ra;
if ( liste[i].C.x > 180. ) liste[i].C.x -= 360.;
if ( liste[i].C.x < -180. ) liste[i].C.x += 360.;
liste[i].C.x *= -3600 * cos(M.ref_dec * DTR);
liste[i].C.x *= -3600 / cos(M.ref_dec * DTR);
liste[i].C.y -= M.ref_dec;
liste[i].C.y *= 3600;
}
......@@ -209,7 +209,7 @@ void getRADEC(char *line,
double ss0, tt0;
int hh0, mm0, dd0, nn0;
char ref1[20], ref2[20];
char ref1[30], ref2[30];
*iref = 0;
*ra = *dec = 0.;
......@@ -245,7 +245,7 @@ void convertXY( double *x, double *y, int iref, double ref_ra, double ref_dec )
// Convert the input values to absolute WCS coordinates
if ( iref == 1 || iref == 3 )
{
*x /= -3600.*cos(ref_dec * DTR);
*x /= -3600./cos(ref_dec * DTR);
*x += ref_ra;
*y /= 3600.;
*y += ref_dec;
......
......@@ -2,12 +2,12 @@
#include<signal.h>
#include<string.h>
#include<math.h>
#include <wcs.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include<lt.h>
/****************************************************************/
/* nom: g_mass */
/* auteur: Jean-Paul Kneib */
......@@ -83,6 +83,7 @@ void g_mass(int imass, int np, double zl, double zs, char *file)
exit(1);
}
// Initialise the mass map(ny,nx)
nx = np;
ny = (int) (nx / (F.xmax - F.xmin) * (F.ymax - F.ymin));
......@@ -127,7 +128,8 @@ void g_mass(int imass, int np, double zl, double zs, char *file)
/* Return a convergence map of size npxnp in pixels
*/
static void computeKmap(int ny, int nx, double dx, double dy, double **map, double zl, double zs)
static void computeKmap(int ny, int nx, double dx, double dy, double **map,
double zl, double zs)
{
const extern struct g_grille G;
const extern struct g_frame F;
......@@ -138,6 +140,8 @@ static void computeKmap(int ny, int nx, double dx, double dy, double **map, doub
int i, j;
long int ilens;
double oldz, dls;
double ra, dec;
int offscl;
// Compute dls for lens[0].z (default for projected kappa)
dls = distcosmo2(lens[0].z, zs);
......@@ -150,6 +154,7 @@ static void computeKmap(int ny, int nx, double dx, double dy, double **map, doub
oldz = zl;
}
#pragma omp parallel for private(i,j,dls,oldz,pi,grad2)
for ( ilens = 0; ilens < G.nlens; ilens ++)
{
if( lens[ilens].z >= zs )
......@@ -161,7 +166,8 @@ static void computeKmap(int ny, int nx, double dx, double dy, double **map, doub
if( zl != 0 && lens[ilens].z != zl )
continue;
printf( "INFO: Compute map for lens %ld/%ld\r", ilens, G.nlens); fflush(stdout);
if( ilens % 4 == 0)
printf( "INFO: Compute map for lens %ld/%ld\r", ilens, G.nlens); fflush(stdout);
// This is for projected kappa
if( lens[ilens].z != oldz )
......
......@@ -110,6 +110,7 @@ void g_shear(int ishear, int np, double z, char *file)
//conv = vol * vol / 4.*distcosmo1(lens[0].z) * dlsds;
//conv = dlsds / 2.;
conv = 1./ dos / 2.;
#pragma omp parallel for private(j,i,grad2)
for (j = 0; j < np; j++)
{
pi.y = j * (ymax - ymin) / (np - 1) + ymin;
......@@ -127,6 +128,7 @@ void g_shear(int ishear, int np, double z, char *file)
//conv = vol * vol / 4.*distcosmo1(lens[0].z) * dlsds;
//conv = dlsds;
conv = 1. / dos;
#pragma omp parallel for private(j,i,grad2)
for (j = 0; j < np; j++)
{
pi.y = j * (ymax - ymin) / (np - 1) + ymin;
......
#include<stdio.h>
#include<math.h>
#include<fonction.h>
#include<constant.h>
#include<dimension.h>
#include<structure.h>
#include<lt.h>
/****************************************************************/
/* nom: g_shear_hp */
/* auteur: Eric Jullo */
/* date: 7/7/2019 */
/* place: San-Francisco */
/****************************************************************/
void g_shear_hp(int ishear, int nside, double z, char *file)
{
const extern struct g_mode M;
const extern struct g_frame F;
const extern struct pot lens[];
register int i, j, ii, jj;
double q, conv;
double dl0s, dos, dlsds;
struct point ps, pi;
struct ellipse ampli;
double **shear;
struct matrix grad2;
double dx, dy, xmin, xmax, ymin, ymax;
long int np;
if (ishear == 1)
{
NPRINTF(stderr, "COMP:shear_map (gamma) in the Image Plane for z_s=%.3lf =>%s\n", z, file);
}
else if (ishear == 2)
{
NPRINTF(stderr, "COMP:shear_map (eps) in the Image Plane for z_s=%.3lf=>%s\n",
z, file);
}
else if (ishear == 3)
{
NPRINTF(stderr, "COMP:shear_map (gamma1) in the Image Plane for z_s=%.3lf=>%s\n",
z, file);
}
else if (ishear == 4)
{
NPRINTF(stderr, "COMP:shear_map (gamma2) in the Image Plane for z_s=%.3lf=>%s\n",
z, file);
}
else if (ishear == -1)
{
NPRINTF(stderr, "COMP:shear_map (gamma) in the Source Plane for z_s=%.3lf =>%s\n",
z, file);
}
else if (ishear == -2)
{
NPRINTF(stderr, "COMP:shear_map (eps) in the Source Plane for z_s=%.3lf =>%s\n",
z, file);
}
// Parse errors
if ( nside == 0 )
{
fprintf(stderr, "ERROR: mass map size set to 0\n");
exit(1);
}
dl0s = distcosmo2( lens[0].z, z );
dos = distcosmo1( z );
dlsds = dl0s / dos;
NPRINTF(stderr, "COMP: Healpix map\n");
// Initialise the mass map(nside)
np = nside2npix(nside);
pixsize = sqrt(4.*PI/12./nside/nside)*RTA;
NPRINTF(stderr, "npix:%ld pixsize:%lf\n", np, pixsize);
shear = (float *) calloc(np, sizeof(float));
// Compute the shear at the positions in the image plane
if (ishear > 0)
{
xmin = F.xmin;
xmax = F.xmax;
ymin = F.ymin;
ymax = F.ymax;
// Compute a shear map
if ( ishear == 1 )
{
for (j = 0; j < np; j++)
{
pi.y = j * (ymax - ymin) / (np - 1) + ymin;
for (i = 0; i < np; i++)
{
pi.x = i * (xmax - xmin) / (np - 1) + xmin;
ampli = e_unmag(&pi, dl0s, dos, z);
shear[j][i] = (ampli.a - ampli.b) / 2.;
}
}
}
// Compute an ellipticity map
else if ( ishear == 2 )
{
for (j = 0; j < np; j++)
{
pi.y = j * (ymax - ymin) / (np - 1) + ymin;
for (i = 0; i < np; i++)
{
pi.x = i * (xmax - xmin) / (np - 1) + xmin;
ampli = e_unmag(&pi, dl0s, dos, z);
q = ampli.a / ampli.b;
shear[j][i] = fabs((q * q - 1.) / (q * q + 1.));
}
}
}
// Compute the gamma1 shear component
else if ( ishear == 3 )
{
//conv = vol * vol / 4.*distcosmo1(lens[0].z) * dlsds;
//conv = dlsds / 2.;
conv = 1./ dos / 2.;
#pragma omp parallel for private(j,i,grad2)
for (j = 0; j < np; j++)
{
pi.y = j * (ymax - ymin) / (np - 1) + ymin;
for (i = 0; i < np; i++)
{
pi.x = i * (xmax - xmin) / (np - 1) + xmin;
grad2 = e_grad2(&pi, dl0s, z);
shear[j][i] = conv * (grad2.a - grad2.c);
}
}
}
// Compute the gamma2 shear component
else if ( ishear == 4 )
{
//conv = vol * vol / 4.*distcosmo1(lens[0].z) * dlsds;
//conv = dlsds;
conv = 1. / dos;
#pragma omp parallel for private(j,i,grad2)
for (i = 0; i < np; i++)
{
pi.x = i * (xmax - xmin) / (np - 1) + xmin;
grad2 = e_grad2(&pi, dl0s, z);
shear[j][i] = conv * grad2.b;
}
}
}
/* else
// Compute the shear at the positions in the source plane
// (what will be felt by each pixel of the source plane)
{
dx = (F.xmax - F.xmin) / 6.;
dy = (F.ymax - F.ymin) / 6.;
xmin = F.xmin + dx; //define a smaller window in the Source plane
xmax = F.xmax - dx;
ymin = F.ymin + dy;
ymax = F.ymax - dy;
dx = (xmax - xmin) / (np - 1);
dy = (ymax - ymin) / (np - 1);
for (j = 0; j < (np*1.8); j++)
{
pi.y = j * (F.ymax - F.ymin) / (np * 1.8 - 1) + F.ymin;
for (i = 0; i < (np*1.8); i++)
{
pi.x = i * (F.xmax - F.xmin) / (np * 1.8 - 1) + F.xmin;
ampli = e_unmag(&pi, dl0s, dos, z);
e_dpl(&pi, dlsds, z, &ps);
ii = (int) (0.5 + (ps.x - xmin) / dx);
jj = (int) (0.5 + (ps.y - ymin) / dy);
if ((ii >= 0) && (ii < np) && (jj >= 0) && (jj < np))
{
if ( ishear == -1)
shear[jj][ii] += (ampli.a - ampli.b) / 2.;
else if ( ishear == -2)
{
q = ampli.a / ampli.b;
shear[jj][ii] += fabs((q * q - 1.) / (q * q + 1.));
}
nsh[jj][ii]++;
}
}
}
for (ii = 0; ii < np; ii++)
for (jj = 0; jj < np; jj++)
if (nsh[jj][ii] > 0)
shear[jj][ii] /= nsh[jj][ii];
}*/
// Write the mass maps
char coordsys[2];
strcpy(coordsys, "E");
remove(file);
write_healpix_map(shear, nside, file, 1, coordsys);
free(shear);
}
......@@ -83,4 +83,6 @@ void r_frame(FILE *IN, FILE *OUT)
if (F.rmax == 0.)
F.rmax = 30.;
F.wcsinfo->nxpix = (int)(2.*F.rmax);
F.wcsinfo->nypix = (int)(2.*F.rmax);
}
......@@ -20,6 +20,7 @@ static void scanInverse(char * third);
void r_runmode(FILE *IN, FILE *OUT)
{
extern struct g_mode M;
extern struct g_frame F;
char second[20], third[FILENAME_SIZE+10];
char ref1[20], ref2[20];
......@@ -263,6 +264,12 @@ void r_runmode(FILE *IN, FILE *OUT)
fprintf(OUT, "\t%s\n", second);
// Set the reference coordinate in the global wcs info of Lenstool
if (M.iref == 3)
wcsreset(F.wcsinfo, F.wcsinfo->crpix[0], F.wcsinfo->crpix[1],
M.ref_ra, M.ref_dec, F.wcsinfo->cdelt[0], F.wcsinfo->cdelt[1], 0.,
F.wcsinfo->cd);
}
static void scanInverse(char * third)
......
......@@ -153,6 +153,11 @@ void set_default()
F.ymin = -20.;
F.ymax = 20.;
F.rmax = 25.;
char ctype1[9], ctype2[9];
strcpy(ctype1, "RA---TAN");
strcpy(ctype2, "DEC--TAN");
F.wcsinfo = wcskinit((int)(2*F.rmax),(int)(2*F.rmax),ctype1,ctype2,
0.,0.,0.,0.,NULL,-1./3600.,1./3600.,0.,2000,0.); // Dummy size but 1arcsec/pixel
/* image plane */
imFrame.xmin = -20.;
......
......@@ -19,10 +19,10 @@ void tirer(struct galaxie *source)
{
const extern struct pot lens[];