diff --git a/kernel/eplp.cl b/kernel/eplp.cl deleted file mode 100644 index 54ba831..0000000 --- a/kernel/eplp.cl +++ /dev/null @@ -1,64 +0,0 @@ -// An elliptical power law potential model \Phi(r,theta) = re ( G(r,\theta)/re )^alpha G(r,\theta) = sqrt( x^2 + q^2 y^2) -// see Leier & Metcalf 2015 -OBJECT(eplp) = LENS; - -PARAMS(eplp) = { - { "x" }, - { "y" }, - { "re" }, - { "alpha" }, - { "q" }, - { "pa", true } -}; - -struct eplp -{ - float2 x; - mat22 m; - mat22 w; - - float q2; - float e; - float re; - float alpha; -}; - -static float2 eplp(constant struct eplp* eplp, float2 x) -{ - float2 dy,y; - float r; - - dy = mv22(eplp->m, x - eplp->x); - - r = length(dy); - - float eta = sqrt(dy.x*dy.x + eplp->q2*dy.y*dy.y); - - float a_iso = eplp->re*powr(eta/eplp->re , eplp->alpha) / eta; - - y.x = a_iso * dy.x ; - y.y = a_iso * eplp->q2 * dy.y; - - return mv22(eplp->w, y); -} - -static void set_eplp(global struct eplp* eplp, float x1, float x2, float re, float alpha, float q, float pa) -{ - float c = cos(pa*DEG2RAD); - float s = sin(pa*DEG2RAD); - - // lens position - eplp->x = (float2)(x1, x2); - - // rotation matrix - eplp->m = (mat22)(c, s, -s, c); - - // inverse rotation matrix - eplp->w = (mat22)(c, -s, s, c); - - eplp->q2 = q*q; - eplp->e = 1 - q*q; - eplp->re = re; - eplp->alpha = alpha; - -} diff --git a/kernel/eplp_plus_shear.cl b/kernel/eplp_plus_shear.cl deleted file mode 100644 index 69df1fd..0000000 --- a/kernel/eplp_plus_shear.cl +++ /dev/null @@ -1,71 +0,0 @@ -// An elliptical power law potential model plus shear -// \Phi(r,theta) = re ( G(r,\theta)/re )^alpha G(r,\theta) = sqrt( x^2 + q^2 y^2) -// see Leier & Metcalf 2015 -OBJECT(eplp_plus_shear) = LENS; - -PARAMS(eplp_plus_shear) = { - { "x" }, - { "y" }, - { "re" }, - { "alpha" }, - { "q" }, - { "pa", true }, - { "g1" }, - { "g2" } - -}; - -struct eplp_plus_shear -{ - float2 x; - mat22 m; - mat22 w; - - float q2; - float e; - float re; - float alpha; - - mat22 g; -}; - -static float2 eplp_plus_shear(constant struct eplp_plus_shear* eplp_plus_shear, float2 x) -{ - float2 y; - - float2 dx = x - eplp_plus_shear->x; - float2 dy = mv22(eplp_plus_shear->m,dx); - - float eta = sqrt(dy.x*dy.x + eplp_plus_shear->q2*dy.y*dy.y); - - float a_iso = eplp_plus_shear->re*powr(eta/eplp_plus_shear->re , eplp_plus_shear->alpha) / eta; - - y.x = a_iso * dy.x ; - y.y = a_iso * eplp_plus_shear->q2 * dy.y; - - y = mv22(eplp_plus_shear->w,y); - - return y + mv22(eplp_plus_shear->g, dx); -} - -static void set_eplp_plus_shear(global struct eplp_plus_shear* eplp_plus_shear, float x1, float x2, float re, float alpha, float q, float pa, float g1, float g2) -{ - float c = cos(pa*DEG2RAD); - float s = sin(pa*DEG2RAD); - - // lens position - eplp_plus_shear->x = (float2)(x1, x2); - - // rotation matrix - eplp_plus_shear->m = (mat22)(c, s, -s, c); - - // inverse rotation matrix - eplp_plus_shear->w = (mat22)(c, -s, s, c); - - eplp_plus_shear->q2 = q*q; - eplp_plus_shear->e = 1 - q*q; - eplp_plus_shear->re = re; - eplp_plus_shear->alpha = alpha; - - eplp_plus_shear->g = (mat22)(g1,g2,g2,-g1); -} diff --git a/kernel/gauss.cl b/kernel/gauss.cl deleted file mode 100644 index 78ed7f0..0000000 --- a/kernel/gauss.cl +++ /dev/null @@ -1,40 +0,0 @@ -OBJECT(gauss) = SOURCE; - -PARAMS(gauss) = { - { "x" }, - { "y" }, - { "sigma" }, - { "mag" }, - { "q" }, - { "pa", true }, -}; - -struct gauss -{ - float2 x; // source position - mat22 t; // coordinate transformation matrix - float s2; // variance - float norm; // normalisation -}; - -static float gauss(constant struct gauss* data, float2 x) -{ - // Gaussian profile for centered and rotated coordinate system - float2 y = mv22(data->t, x - data->x); - return data->norm*exp(-0.5f*dot(y, y)/data->s2); -} - -static void set_gauss(global struct gauss* data, float x1, float x2, float sigma, float mag, float q, float pa) -{ - float c = cos(pa*DEG2RAD); - float s = sin(pa*DEG2RAD); - - // source position - data->x = (float2)(x1, x2); - - // transformation matrix: rotate and scale - data->t = (mat22)(q*c, q*s, -s, c); - - data->s2 = sigma*sigma; - data->norm = exp(-0.4f*mag*LOG_10)*0.5f/PI/data->s2/q; -} diff --git a/kernel/object.cl b/kernel/object.cl index 6c57637..d56804e 100644 --- a/kernel/object.cl +++ b/kernel/object.cl @@ -12,15 +12,3 @@ struct __attribute__ ((aligned(4))) param char name[32]; int wrap; }; - -// macro to specify type of object -#define OBJECT(x) constant int object_##x - -// macro to simplify parameter definition for object -#define PARAMS(x) constant struct param parlst_##x[] - -// macro to get number of parameters for object -#define NPARAMS(x) (sizeof(parlst_##x)/sizeof(struct param)) - -// macro to access specific parameter definition for object -#define PARAM(x, n) parlst_##x[n] diff --git a/kernel/point_mass.cl b/kernel/point_mass.cl deleted file mode 100644 index 7b27361..0000000 --- a/kernel/point_mass.cl +++ /dev/null @@ -1,30 +0,0 @@ -// point mass lens - -OBJECT(point_mass) = LENS; - -PARAMS(point_mass) = { - { "x" }, - { "y" }, - { "r" } -}; - -struct point_mass -{ - float2 x; // lens position - float r2; // Einstein radius squared -}; - -static float2 point_mass(constant struct point_mass* data, float2 x) -{ - // point mass deflection - return data->r2*normalize(x - data->x); -} - -static void set_point_mass(global struct point_mass* data, float x1, float x2, float r) -{ - // lens position - data->x = (float2)(x1, x2); - - // Einstein radius squared - data->r2 = r*r; -} diff --git a/kernel/sersic.cl b/kernel/sersic.cl deleted file mode 100644 index ea9b76b..0000000 --- a/kernel/sersic.cl +++ /dev/null @@ -1,44 +0,0 @@ -OBJECT(sersic) = SOURCE; - -PARAMS(sersic) = { - { "x" }, - { "y" }, - { "r" }, - { "mag" }, - { "n" }, - { "q" }, - { "pa", true }, -}; - -struct sersic -{ - float2 x; // source position - mat22 t; // coordinate transformation matrix - float log0; // profile constants - float log1; - float m; -}; - -static float sersic(constant struct sersic* data, float2 x) -{ - float2 y = mv22(data->t, x - data->x); - return exp(data->log0 - exp(data->log1 + data->m*log(dot(y, y)))); -} - -static void set_sersic(global struct sersic* data, float x1, float x2, float r, float mag, float n, float q, float pa) -{ - float b = 1.9992f*n - 0.3271f; // approximation valid for 0.5 < n < 8 - - float c = cos(pa*DEG2RAD); - float s = sin(pa*DEG2RAD); - - // source position - data->x = (float2)(x1, x2); - - // transformation matrix: rotate and scale - data->t = (mat22)(q*c, q*s, -s, c); - - data->log0 = -0.4f*mag*LOG_10 + 2*n*log(b) - LOG_PI - 2*log(r) - log(tgamma(2*n+1)); - data->log1 = log(b) - (0.5f*log(q) + log(r))/n; - data->m = 0.5f/n; -} diff --git a/kernel/sis_plus_shear.cl b/kernel/sis_plus_shear.cl deleted file mode 100644 index 6dcd0ae..0000000 --- a/kernel/sis_plus_shear.cl +++ /dev/null @@ -1,39 +0,0 @@ -// singular isothermal sphere plus shear - -OBJECT(sis_plus_shear) = LENS; - -PARAMS(sis_plus_shear) = { - { "x" }, - { "y" }, - { "r" }, - { "g1" }, - { "g2" } -}; - -struct sis_plus_shear -{ - float2 x; // lens position - mat22 g; // shear matrix - float r; // Einstein radius -}; - -static float2 sis_plus_shear(constant struct sis_plus_shear* data, float2 x) -{ - // move to central coordinates - x -= data->x; - - // SIS deflection plus external shear - return data->r*normalize(x) + mv22(data->g, x); -} - -static void set_sis_plus_shear(global struct sis_plus_shear* data, float x1, float x2, float r, float g1, float g2) -{ - // lens position - data->x = (float2)(x1, x2); - - // Einstein radius - data->r = r; - - // shear matrix - data->g = (mat22)(g1, g2, g2, -g1); -} diff --git a/kernel/sky.cl b/kernel/sky.cl deleted file mode 100644 index 477269b..0000000 --- a/kernel/sky.cl +++ /dev/null @@ -1,24 +0,0 @@ -OBJECT(sky) = FOREGROUND; - -PARAMS(sky) = { - { "bg" }, - { "dx" }, - { "dy" } -}; - -struct sky -{ - float bg; - float2 grad; -}; - -static float sky(constant struct sky* sky, float2 x) -{ - return sky->bg + dot(sky->grad, x - (float2)(1, 1)); -} - -static void set_sky(global struct sky* sky, float bg, float dx, float dy) -{ - sky->bg = bg; - sky->grad = (float2)(dx, dy); -} diff --git a/kernel/devauc.cl b/objects/devauc.cl similarity index 50% rename from kernel/devauc.cl rename to objects/devauc.cl index b934aaa..c384c7e 100644 --- a/kernel/devauc.cl +++ b/objects/devauc.cl @@ -1,11 +1,12 @@ -// De Vaucouleurs b param +// de Vaucouleurs b param #define DEVAUC_B 7.6692494425008039044f // (DEVAUC_B^8)/(8!) #define DEVAUC_C 296.826303766893f -OBJECT(devauc) = SOURCE; +type = SOURCE; -PARAMS(devauc) = { +params +{ { "x" }, { "y" }, { "r" }, @@ -14,7 +15,7 @@ PARAMS(devauc) = { { "pa", true }, }; -struct devauc +data { float2 x; // source position mat22 t; // coordinate transformation matric @@ -22,26 +23,26 @@ struct devauc float norm; // normalisation }; -static float devauc(constant struct devauc* data, float2 x) +static float brightness(constant data* this, float2 x) { - // DeVaucouleur's profile for centered and rotated coordinate system - return data->norm*exp(-DEVAUC_B*sqrt(sqrt(length(mv22(data->t, x - data->x))/data->rs))); + // de Vaucouleurs profile for centered and rotated coordinate system + return this->norm*exp(-DEVAUC_B*sqrt(sqrt(length(mv22(this->t, x - this->x))/this->rs))); } -static void set_devauc(global struct devauc* data, float x1, float x2, float r, float mag, float q, float pa) +static void set(global data* this, float x, float y, float r, float mag, float q, float pa) { float c = cos(pa*DEG2RAD); float s = sin(pa*DEG2RAD); // source position - data->x = (float2)(x1, x2); + this->x = (float2)(x, y); // transformation matrix: rotate and scale - data->t = (mat22)(q*c, q*s, -s, c); + this->t = (mat22)(q*c, q*s, -s, c); // scale length - data->rs = r; + this->rs = r; // normalisation to total luminosity - data->norm = exp(-0.4f*mag*LOG_10)/PI/r/r/q*DEVAUC_C; + this->norm = exp(-0.4f*mag*LOG_10)/PI/r/r/q*DEVAUC_C; } diff --git a/objects/eplp.cl b/objects/eplp.cl new file mode 100644 index 0000000..9b8af05 --- /dev/null +++ b/objects/eplp.cl @@ -0,0 +1,64 @@ +// An elliptical power law potential model \Phi(r,theta) = re ( G(r,\theta)/re )^alpha G(r,\theta) = sqrt( x^2 + q^2 y^2) +// see Leier & Metcalf 2015 +type = LENS; + +params +{ + { "x" }, + { "y" }, + { "re" }, + { "alpha" }, + { "q" }, + { "pa", true } +}; + +data +{ + float2 x; + mat22 m; + mat22 w; + + float q2; + float e; + float re; + float alpha; +}; + +static float2 deflection(constant data* this, float2 x) +{ + float2 dy,y; + float r; + + dy = mv22(this->m, x - this->x); + + r = length(dy); + + float eta = sqrt(dy.x*dy.x + this->q2*dy.y*dy.y); + + float a_iso = this->re*powr(eta/this->re , this->alpha) / eta; + + y.x = a_iso * dy.x ; + y.y = a_iso * this->q2 * dy.y; + + return mv22(this->w, y); +} + +static void set(global data* this, float x, float y, float re, float alpha, float q, float pa) +{ + float c = cos(pa*DEG2RAD); + float s = sin(pa*DEG2RAD); + + // lens position + this->x = (float2)(x, y); + + // rotation matrix + this->m = (mat22)(c, s, -s, c); + + // inverse rotation matrix + this->w = (mat22)(c, -s, s, c); + + this->q2 = q*q; + this->e = 1 - q*q; + this->re = re; + this->alpha = alpha; +} diff --git a/objects/eplp_plus_shear.cl b/objects/eplp_plus_shear.cl new file mode 100644 index 0000000..08e39f5 --- /dev/null +++ b/objects/eplp_plus_shear.cl @@ -0,0 +1,71 @@ +// An elliptical power law potential model plus shear +// \Phi(r,theta) = re ( G(r,\theta)/re )^alpha G(r,\theta) = sqrt( x^2 + q^2 y^2) +// see Leier & Metcalf 2015 +type = LENS; + +params +{ + { "x" }, + { "y" }, + { "re" }, + { "alpha" }, + { "q" }, + { "pa", true }, + { "g1" }, + { "g2" } +}; + +data +{ + float2 x; + mat22 m; + mat22 w; + + float q2; + float e; + float re; + float alpha; + + mat22 g; +}; + +static float2 deflection(constant data* this, float2 x) +{ + float2 y; + + float2 dx = x - this->x; + float2 dy = mv22(this->m,dx); + + float eta = sqrt(dy.x*dy.x + this->q2*dy.y*dy.y); + + float a_iso = this->re*powr(eta/this->re , this->alpha) / eta; + + y.x = a_iso * dy.x ; + y.y = a_iso * this->q2 * dy.y; + + y = mv22(this->w,y); + + return y + mv22(this->g, dx); +} + +static void set(global data* this, float x, float y, float re, float alpha, float q, float pa, float g1, float g2) +{ + float c = cos(pa*DEG2RAD); + float s = sin(pa*DEG2RAD); + + // lens position + this->x = (float2)(x, y); + + // rotation matrix + this->m = (mat22)(c, s, -s, c); + + // inverse rotation matrix + this->w = (mat22)(c, -s, s, c); + + this->q2 = q*q; + this->e = 1 - q*q; + this->re = re; + this->alpha = alpha; + + this->g = (mat22)(g1,g2,g2,-g1); +} diff --git a/kernel/exponential.cl b/objects/exponential.cl similarity index 52% rename from kernel/exponential.cl rename to objects/exponential.cl index 8f94d26..cc19431 100644 --- a/kernel/exponential.cl +++ b/objects/exponential.cl @@ -1,6 +1,7 @@ -OBJECT(exponential) = SOURCE; +type = SOURCE; -PARAMS(exponential) = { +params +{ { "x" }, { "y" }, { "rs" }, @@ -9,7 +10,7 @@ PARAMS(exponential) = { { "pa", true }, }; -struct exponential +data { float2 x; // source position mat22 t; // coordinate transformation matrix @@ -17,26 +18,26 @@ struct exponential float norm; // normalisation }; -static float exponential(constant struct exponential* data, float2 x) +static float brightness(constant data* this, float2 x) { // exponential profile for centered and rotated coordinate system - return data->norm*exp(-length(mv22(data->t, x - data->x))/data->rs); + return this->norm*exp(-length(mv22(this->t, x - this->x))/this->rs); } -static void set_exponential(global struct exponential* data, float x1, float x2, float rs, float mag, float q, float pa) +static void set(global data* this, float x, float y, float rs, float mag, float q, float pa) { float c = cos(pa*DEG2RAD); float s = sin(pa*DEG2RAD); // source position - data->x = (float2)(x1, x2); + this->x = (float2)(x, y); // transformation matrix: rotate and scale - data->t = (mat22)(q*c, q*s, -s, c); + this->t = (mat22)(q*c, q*s, -s, c); // scale length - data->rs = rs; + this->rs = rs; // normalisation to total luminosity - data->norm = exp(-0.4f*mag*LOG_10)*0.5f/PI/rs/rs/q; + this->norm = exp(-0.4f*mag*LOG_10)*0.5f/PI/rs/rs/q; } diff --git a/objects/gauss.cl b/objects/gauss.cl new file mode 100644 index 0000000..8a89a85 --- /dev/null +++ b/objects/gauss.cl @@ -0,0 +1,41 @@ +type = SOURCE; + +params +{ + { "x" }, + { "y" }, + { "sigma" }, + { "mag" }, + { "q" }, + { "pa", true }, +}; + +data +{ + float2 x; // source position + mat22 t; // coordinate transformation matrix + float s2; // variance + float norm; // normalisation +}; + +static float brightness(constant data* this, float2 x) +{ + // Gaussian profile for centered and rotated coordinate system + float2 y = mv22(this->t, x - this->x); + return this->norm*exp(-0.5f*dot(y, y)/this->s2); +} + +static void set(global data* this, float x, float y, float sigma, float mag, float q, float pa) +{ + float c = cos(pa*DEG2RAD); + float s = sin(pa*DEG2RAD); + + // source position + this->x = (float2)(x, y); + + // transformation matrix: rotate and scale + this->t = (mat22)(q*c, q*s, -s, c); + + this->s2 = sigma*sigma; + this->norm = exp(-0.4f*mag*LOG_10)*0.5f/PI/this->s2/q; +} diff --git a/kernel/nsie.cl b/objects/nsie.cl similarity index 55% rename from kernel/nsie.cl rename to objects/nsie.cl index 97b12bc..3cf64f3 100644 --- a/kernel/nsie.cl +++ b/objects/nsie.cl @@ -1,9 +1,10 @@ // non-singular isothermal ellipsoid // follows Schneider, Kochanek, Wambsganss (2006) -OBJECT(nsie) = LENS; +type = LENS; -PARAMS(nsie) = { +params +{ { "x" }, { "y" }, { "r" }, @@ -12,7 +13,7 @@ PARAMS(nsie) = { { "pa", true } }; -struct nsie +data { float2 x; // lens position mat22 m; // rotation matrix for position angle @@ -25,44 +26,44 @@ struct nsie float d; }; -static float2 nsie(constant struct nsie* data, float2 x) +static float2 deflection(constant data* this, float2 x) { float2 y; float r; // move to central coordinates - x -= data->x; + x -= this->x; // rotate coordinates by position angle - y = mv22(data->m, x); + y = mv22(this->m, x); // NSIE deflection - r = sqrt(data->q2*y.x*y.x + y.y*y.y); - y = data->d*(float2)(atan(y.x*data->e/(data->rc + r)), atanh(y.y*data->e/(data->rc*data->q2 + r))); + r = sqrt(this->q2*y.x*y.x + y.y*y.y); + y = this->d*(float2)(atan(y.x*this->e/(this->rc + r)), atanh(y.y*this->e/(this->rc*this->q2 + r))); // reverse coordinate rotation - return mv22(data->w, y); + return mv22(this->w, y); } -static void set_nsie(global struct nsie* data, float x1, float x2, float r, float rc, float q, float pa) +static void set(global data* this, float x, float y, float r, float rc, float q, float pa) { float c = cos(pa*DEG2RAD); float s = sin(pa*DEG2RAD); // lens position - data->x = (float2)(x1, x2); + this->x = (float2)(x, y); // rotation matrix - data->m = (mat22)(c, s, -s, c); + this->m = (mat22)(c, s, -s, c); // inverse rotation matrix - data->w = (mat22)(c, -s, s, c); + this->w = (mat22)(c, -s, s, c); // core radius - data->rc = rc; + this->rc = rc; // auxiliary quantities - data->q2 = q*q; - data->e = sqrt(1 - q*q); - data->d = r*sqrt(q)/sqrt(1 - q*q); + this->q2 = q*q; + this->e = sqrt(1 - q*q); + this->d = r*sqrt(q)/sqrt(1 - q*q); } diff --git a/kernel/nsis.cl b/objects/nsis.cl similarity index 54% rename from kernel/nsis.cl rename to objects/nsis.cl index ab6f195..e46cda1 100644 --- a/kernel/nsis.cl +++ b/objects/nsis.cl @@ -1,39 +1,40 @@ // non-singular isothermal sphere // follows Schneider, Kochanek, Wambsganss (2006) -OBJECT(nsis) = LENS; +type = LENS; -PARAMS(nsis) = { +params +{ { "x" }, { "y" }, { "r" }, { "rc" } }; -struct nsis +data { float2 x; // lens position float r; // Einstein radius float rc; // core radius }; -static float2 nsis(constant struct nsis* data, float2 x) +static float2 deflection(constant data* this, float2 x) { // move to central coordinates - x -= data->x; + x -= this->x; // NSIS deflection - return data->r/(data->rc + length(x))*x; + return this->r/(this->rc + length(x))*x; } -static void set_nsis(global struct nsis* nsis, float x1, float x2, float r, float rc) +static void set(global data* this, float x, float y, float r, float rc) { // lens position - nsis->x = (float2)(x1, x2); + this->x = (float2)(x, y); // Einstein radius - nsis->r = r; + this->r = r; // core radius - nsis->rc = rc; + this->rc = rc; } diff --git a/objects/point_mass.cl b/objects/point_mass.cl new file mode 100644 index 0000000..b291cf3 --- /dev/null +++ b/objects/point_mass.cl @@ -0,0 +1,32 @@ +// point mass lens + +type = LENS; + +params +{ + { "x" }, + { "y" }, + { "r" } +}; + +data +{ + float2 x; // lens position + float r2; // Einstein radius squared +}; + +static float2 deflection(constant data* this, float2 x) +{ + // point mass deflection + x -= this->x; + return this->r2/dot(x, x)*x; +} + +static void set(global data* this, float x, float y, float r) +{ + // lens position + this->x = (float2)(x, y); + + // Einstein radius squared + this->r2 = r*r; +} diff --git a/objects/sersic.cl b/objects/sersic.cl new file mode 100644 index 0000000..24a5617 --- /dev/null +++ b/objects/sersic.cl @@ -0,0 +1,45 @@ +type = SOURCE; + +params +{ + { "x" }, + { "y" }, + { "r" }, + { "mag" }, + { "n" }, + { "q" }, + { "pa", true } +}; + +data +{ + float2 x; // source position + mat22 t; // coordinate transformation matrix + float log0; // profile constants + float log1; + float m; +}; + +static float brightness(constant data* this, float2 x) +{ + float2 y = mv22(this->t, x - this->x); + return exp(this->log0 - exp(this->log1 + this->m*log(dot(y, y)))); +} + +static void set(global data* this, float x, float y, float r, float mag, float n, float q, float a) +{ + float b = 1.9992f*n - 0.3271f; // approximation valid for 0.5 < n < 8 + + float c = cos(a*DEG2RAD); + float s = sin(a*DEG2RAD); + + // source position + this->x = (float2)(x, y); + + // transformation matrix: rotate and scale + this->t = (mat22)(q*c, q*s, -s, c); + + this->log0 = -0.4f*mag*LOG_10 + 2*n*log(b) - LOG_PI - 2*log(r) - log(tgamma(2*n+1)); + this->log1 = log(b) - (0.5f*log(q) + log(r))/n; + this->m = 0.5f/n; +} diff --git a/kernel/sie.cl b/objects/sie.cl similarity index 56% rename from kernel/sie.cl rename to objects/sie.cl index 140791d..c97b8fa 100644 --- a/kernel/sie.cl +++ b/objects/sie.cl @@ -1,9 +1,10 @@ // singular isothermal ellipsoid // follows Schneider, Kochanek, Wambsganss (2006) -OBJECT(sie) = LENS; +type = LENS; -PARAMS(sie) = { +params +{ { "x" }, { "y" }, { "r" }, @@ -11,7 +12,7 @@ PARAMS(sie) = { { "pa", true } }; -struct sie +data { float2 x; // lens position mat22 m; // rotation matrix for position angle @@ -23,41 +24,41 @@ struct sie float d; }; -static float2 sie(constant struct sie* data, float2 x) +static float2 deflection(constant data* this, float2 x) { float2 y; float r; // move to central coordinates - x -= data->x; + x -= this->x; // rotate coordinates by position angle - y = mv22(data->m, x); + y = mv22(this->m, x); // SIE deflection - r = data->e/sqrt(data->q2*y.x*y.x + y.y*y.y); - y = data->d*(float2)(atan(y.x*r), atanh(y.y*r)); + r = this->e/sqrt(this->q2*y.x*y.x + y.y*y.y); + y = this->d*(float2)(atan(y.x*r), atanh(y.y*r)); // reverse coordinate rotation - return mv22(data->w, y); + return mv22(this->w, y); } -static void set_sie(global struct sie* sie, float x1, float x2, float r, float q, float pa) +static void set(global data* this, float x, float y, float r, float q, float pa) { float c = cos(pa*DEG2RAD); float s = sin(pa*DEG2RAD); // lens position - sie->x = (float2)(x1, x2); + this->x = (float2)(x, y); // rotation matrix - sie->m = (mat22)(c, s, -s, c); + this->m = (mat22)(c, s, -s, c); // inverse rotation matrix - sie->w = (mat22)(c, -s, s, c); + this->w = (mat22)(c, -s, s, c); // auxiliary quantities - sie->q2 = q*q; - sie->e = sqrt(1 - q*q); - sie->d = r*sqrt(q)/sqrt(1 - q*q); + this->q2 = q*q; + this->e = sqrt(1 - q*q); + this->d = r*sqrt(q)/sqrt(1 - q*q); } diff --git a/kernel/sie_plus_shear.cl b/objects/sie_plus_shear.cl similarity index 53% rename from kernel/sie_plus_shear.cl rename to objects/sie_plus_shear.cl index 187c5db..0dc0b59 100644 --- a/kernel/sie_plus_shear.cl +++ b/objects/sie_plus_shear.cl @@ -1,9 +1,10 @@ // singular isothermal ellipsoid plus shear // follows Schneider, Kochanek, Wambsganss (2006) -OBJECT(sie_plus_shear) = LENS; +type = LENS; -PARAMS(sie_plus_shear) = { +params +{ { "x" }, { "y" }, { "r" }, @@ -13,7 +14,7 @@ PARAMS(sie_plus_shear) = { { "g2" } }; -struct sie_plus_shear +data { float2 x; // lens position mat22 m; // rotation matrix for position angle @@ -26,44 +27,44 @@ struct sie_plus_shear float d; }; -static float2 sie_plus_shear(constant struct sie_plus_shear* data, float2 x) +static float2 deflection(constant data* this, float2 x) { float2 y; float r; // move to central coordinates - x -= data->x; + x -= this->x; // rotate coordinates by position angle - y = mv22(data->m, x); + y = mv22(this->m, x); // SIE deflection - r = data->e/sqrt(data->q2*y.x*y.x + y.y*y.y); - y = data->d*(float2)(atan(y.x*r), atanh(y.y*r)); + r = this->e/sqrt(this->q2*y.x*y.x + y.y*y.y); + y = this->d*(float2)(atan(y.x*r), atanh(y.y*r)); // reverse coordinate rotation, apply shear - return mv22(data->w, y) + mv22(data->g, x); + return mv22(this->w, y) + mv22(this->g, x); } -static void set_sie_plus_shear(global struct sie_plus_shear* data, float x1, float x2, float r, float q, float pa, float g1, float g2) +static void set(global data* this, float x, float y, float r, float q, float pa, float g1, float g2) { float c = cos(pa*DEG2RAD); float s = sin(pa*DEG2RAD); // lens position - data->x = (float2)(x1, x2); + this->x = (float2)(x, y); // rotation matrix - data->m = (mat22)(c, s, -s, c); + this->m = (mat22)(c, s, -s, c); // inverse rotation matrix - data->w = (mat22)(c, -s, s, c); + this->w = (mat22)(c, -s, s, c); // shear matrix - data->g = (mat22)(g1, g2, g2, -g1); + this->g = (mat22)(g1, g2, g2, -g1); // auxiliary quantities - data->q2 = q*q; - data->e = sqrt(1 - q*q); - data->d = r*sqrt(q)/sqrt(1 - q*q); + this->q2 = q*q; + this->e = sqrt(1 - q*q); + this->d = r*sqrt(q)/sqrt(1 - q*q); } diff --git a/kernel/sis.cl b/objects/sis.cl similarity index 50% rename from kernel/sis.cl rename to objects/sis.cl index 3a35a37..8b9f165 100644 --- a/kernel/sis.cl +++ b/objects/sis.cl @@ -1,31 +1,32 @@ // singular isothermal sphere // follows Schneider, Kochanek, Wambsganss (2006) -OBJECT(sis) = LENS; +type = LENS; -PARAMS(sis) = { +params +{ { "x" }, { "y" }, { "r" } }; -struct sis +data { float2 x; // lens position float r; // Einstein radius }; -static float2 sis(constant struct sis* data, float2 x) +static float2 deflection(constant data* this, float2 x) { // SIS deflection - return data->r*normalize(x - data->x); + return this->r*normalize(x - this->x); } -static void set_sis(global struct sis* data, float x1, float x2, float r) +static void set(global data* this, float x, float y, float r) { // lens position - data->x = (float2)(x1, x2); + this->x = (float2)(x, y); // Einstein radius - data->r = r; + this->r = r; } diff --git a/objects/sis_plus_shear.cl b/objects/sis_plus_shear.cl new file mode 100644 index 0000000..d79f0ee --- /dev/null +++ b/objects/sis_plus_shear.cl @@ -0,0 +1,40 @@ +// singular isothermal sphere plus shear + +type = LENS; + +params +{ + { "x" }, + { "y" }, + { "r" }, + { "g1" }, + { "g2" } +}; + +data +{ + float2 x; // lens position + mat22 g; // shear matrix + float r; // Einstein radius +}; + +static float2 deflection(constant data* this, float2 x) +{ + // move to central coordinates + x -= this->x; + + // SIS deflection plus external shear + return this->r*normalize(x) + mv22(this->g, x); +} + +static void set(global data* this, float x, float y, float r, float g1, float g2) +{ + // lens position + this->x = (float2)(x, y); + + // Einstein radius + this->r = r; + + // shear matrix + this->g = (mat22)(g1, g2, g2, -g1); +} diff --git a/objects/sky.cl b/objects/sky.cl new file mode 100644 index 0000000..a0ec63a --- /dev/null +++ b/objects/sky.cl @@ -0,0 +1,25 @@ +type = FOREGROUND; + +params +{ + { "bg" }, + { "dx" }, + { "dy" } +}; + +data +{ + float bg; + float2 grad; +}; + +static float foreground(constant data* this, float2 x) +{ + return this->bg + dot(this->grad, x - (float2)(1, 1)); +} + +static void set(global data* this, float bg, float dx, float dy) +{ + this->bg = bg; + this->grad = (float2)(dx, dy); +} diff --git a/src/kernel.c b/src/kernel.c index d62926f..4578d23 100644 --- a/src/kernel.c +++ b/src/kernel.c @@ -8,14 +8,15 @@ #include "kernel.h" #include "log.h" -// parts of kernel files +// parts of kernel and object files static const char KERNEL_DIR[] = "kernel/"; -static const char KERNEL_EXT[] = ".cl"; +static const char OBJECT_DIR[] = "objects/"; +static const char OPENCL_EXT[] = ".cl"; // marker for individual files in kernel code static const char FILEHEAD[] = "//----------------------------------------------------------------------------\n" - "// %s\n" + "// %s%s\n" "//----------------------------------------------------------------------------\n" "\n" ; @@ -38,24 +39,24 @@ static const size_t NMAINKERNS = sizeof(MAINKERNS)/sizeof(MAINKERNS[0]); // kernel to get meta-data for object static const char METAKERN[] = - "kernel void meta_(global int* type, global ulong* size, global ulong* npars)\n" + "kernel void meta_(global int* type, global ulong* size, global ulong* npar)\n" "{\n" - " *type = object_;\n" - " *size = sizeof(struct );\n" - " *npars = NPARAMS();\n" + " *type = type_;\n" + " *size = sizeof(struct data_);\n" + " *npar = sizeof(parlst_)/sizeof(struct param);\n" "}\n" ; // kernel to get parameters for object static const char PARSKERN[] = - "kernel void params_(global struct param* params)\n" + "kernel void params_(global struct param* par)\n" "{\n" - " for(size_t i = 0; i < NPARAMS(); ++i)\n" + " for(size_t i = 0; i < sizeof(parlst_)/sizeof(struct param); ++i)\n" " {\n" - " for(size_t j = 0; j < sizeof(params[i].name); ++j)\n" - " params[i].name[j] = PARAM(, i).name[j];\n" + " for(size_t j = 0; j < sizeof(par[i].name); ++j)\n" + " par[i].name[j] = parlst_[i].name[j];\n" " \n" - " params[i].wrap = PARAM(, i).wrap;\n" + " par[i].wrap = parlst_[i].wrap;\n" " }\n" "}\n" ; @@ -64,40 +65,43 @@ static const char PARSKERN[] = static const char COMPHEAD[] = "static float compute(constant char* data, float2 x)\n" "{\n" - " // initial ray position\n" + " // ray position\n" " float2 y = x;\n" " \n" - " // initial deflection is zero\n" - " float2 a = 0;\n" - " \n" " // initial surface brightness is zero\n" " float f = 0;\n" ; static const char COMPLHED[] = " \n" - " // calculate deflection\n" + " // lens plane\n" + " {\n" + " // initial deflection is zero\n" + " float2 a = 0;\n" + " \n" + " // calculate deflection\n" ; static const char COMPLENS[] = - " a += %s((constant void*)(data + %zu), y);\n" + " a += deflection_%s((constant void*)(data + %zu), y);\n" ; static const char COMPDEFL[] = - " \n" - " // apply deflection to ray, if finite\n" - " y -= dot(a,a) < HUGE_VALF ? a : (float2)(1E10f, 1E10f);\n" + " \n" + " // apply deflection to ray, if finite\n" + " y -= dot(a,a) < HUGE_VALF ? a : (float2)(1E10f, 1E10f);\n" + " }\n" ; static const char COMPSHED[] = " \n" " // calculate surface brightness\n" ; static const char COMPSRCE[] = - " f += %s((constant void*)(data + %zu), y);\n" + " f += brightness_%s((constant void*)(data + %zu), y);\n" ; static const char COMPFHED[] = " \n" " // add foreground\n" ; static const char COMPFGND[] = - " f += %s((constant void*)(data + %zu), x);\n" + " f += foreground_%s((constant void*)(data + %zu), x);\n" ; static const char COMPFOOT[] = " \n" @@ -118,6 +122,28 @@ static const char SETPFOOT[] = "}\n" ; +// object kernel +static const char OBJHEAD[] = + "#define type constant int type_%s\n" + "#define params constant struct param parlst_%s[] = \n" + "#define data struct data_%s\n" + "#define deflection deflection_%s\n" + "#define brightness brightness_%s\n" + "#define foreground foreground_%s\n" + "#define set set_%s\n" + "\n" +; +static const char OBJFOOT[] = + "\n" + "#undef type\n" + "#undef params\n" + "#undef data\n" + "#undef deflection\n" + "#undef brightness\n" + "#undef foreground\n" + "#undef set\n" +; + // replace substring, used to fill in object names in kernels static const char* str_replace(const char* str, const char* search, const char* replace) { @@ -251,7 +277,7 @@ static const char* compute_kernel(size_t nobjs, object objs[]) out = buf; // write file header - wri = sprintf(out, FILEHEAD, "compute"); + wri = sprintf(out, FILEHEAD, "", "compute"); if(wri < 0) errori(NULL); out += wri; @@ -382,7 +408,7 @@ static const char* set_params_kernel(size_t nobjs, object objs[]) out = buf; // write file header - wri = sprintf(out, FILEHEAD, "set_params"); + wri = sprintf(out, FILEHEAD, "", "set_params"); if(wri < 0) errori(NULL); out += wri; @@ -456,8 +482,8 @@ static const char* load_kernel(const char* name) int wri; // construct filename for kernel - filename = malloc(strlen(LENSED_PATH) + strlen(KERNEL_DIR) + strlen(name) + strlen(KERNEL_EXT) + 1); - sprintf(filename, "%s%s%s%s", LENSED_PATH, KERNEL_DIR, name, KERNEL_EXT); + filename = malloc(strlen(LENSED_PATH) + strlen(KERNEL_DIR) + strlen(name) + strlen(OPENCL_EXT) + 1); + sprintf(filename, "%s%s%s%s", LENSED_PATH, KERNEL_DIR, name, OPENCL_EXT); // try to read file f = fopen(filename, "r"); @@ -472,10 +498,12 @@ static const char* load_kernel(const char* name) file_size = ftell(f); // calculate size of buffer - buf_size = file_size + sizeof(FILEHEAD) + strlen(name) + sizeof(FILEFOOT); + buf_size = file_size + + snprintf(NULL, 0, FILEHEAD, KERNEL_DIR, name) + + snprintf(NULL, 0, FILEFOOT); // try to allocate buffer - buf = malloc(buf_size); + buf = malloc(buf_size + 1); if(!buf) errori("kernel %s", name); @@ -483,7 +511,7 @@ static const char* load_kernel(const char* name) out = buf; // write file header - wri = sprintf(out, FILEHEAD, name); + wri = sprintf(out, FILEHEAD, KERNEL_DIR, name); if(wri < 0) errori("kernel %s", name); out += wri; @@ -507,6 +535,91 @@ static const char* load_kernel(const char* name) return buf; } +static const char* load_object(const char* name) +{ + // file for object + char* filename; + FILE* f; + long file_size; + + // buffer for object code + char* buf; + size_t buf_size; + + // current output position + char* out; + + // number of characters added + int wri; + + // construct filename for object + filename = malloc(strlen(LENSED_PATH) + strlen(OBJECT_DIR) + strlen(name) + strlen(OPENCL_EXT) + 1); + sprintf(filename, "%s%s%s%s", LENSED_PATH, OBJECT_DIR, name, OPENCL_EXT); + + // try to read file + f = fopen(filename, "r"); + if(!f) + { + verbose("file not found: %s", filename); + error("could not load object \"%s\"", name); + } + + // go to end of file and get its size + fseek(f, 0, SEEK_END); + file_size = ftell(f); + + // calculate size of buffer + buf_size = file_size + + snprintf(NULL, 0, FILEHEAD, OBJECT_DIR, name) + + snprintf(NULL, 0, OBJHEAD, name, name, name, name, name, name, name) + + snprintf(NULL, 0, OBJFOOT) + + snprintf(NULL, 0, FILEFOOT); + + // try to allocate buffer + buf = malloc(buf_size + 1); + if(!buf) + errori("object %s", name); + + // start at beginning + out = buf; + + // write file header + wri = sprintf(out, FILEHEAD, OBJECT_DIR, name); + if(wri < 0) + errori("object %s", name); + out += wri; + + // write object header + wri = sprintf(out, OBJHEAD, name, name, name, name, name, name, name); + if(wri < 0) + errori("object %s", name); + out += wri; + + // write file + fseek(f, 0, SEEK_SET); + fread(out, 1, file_size, f); + out += file_size; + + // write object footer + wri = sprintf(out, OBJFOOT); + if(wri < 0) + errori("object %s", name); + out += wri; + + // write file footer + wri = sprintf(out, FILEFOOT); + if(wri < 0) + errori("object %s", name); + out += wri; + + // clean up + fclose(f); + free(filename); + + // done + return buf; +} + void object_program(const char* name, size_t* nkernels, const char*** kernels) { // create kernel array with space for object and system kernels @@ -520,7 +633,7 @@ void object_program(const char* name, size_t* nkernels, const char*** kernels) *(k++) = load_kernel(INITKERNS[i]); // load kernel for object - *(k++) = load_kernel(name); + *(k++) = load_object(name); // add kernels for object meta-data and parameters *(k++) = str_replace(METAKERN, "", name); @@ -554,7 +667,7 @@ void main_program(size_t nobjs, object objs[], size_t* nkernels, const char*** k // load kernels for objects for(size_t i = 0; i < nuniq; ++i) - *(k++) = load_kernel(uniq[i]); + *(k++) = load_object(uniq[i]); // load compute kernel *(k++) = compute_kernel(nobjs, objs); diff --git a/tests/lens/point_mass.fits b/tests/lens/point_mass.fits index 7fd1d98..9b9bcbb 100644 Binary files a/tests/lens/point_mass.fits and b/tests/lens/point_mass.fits differ