pro shapelets_shear, structure, gammas, $
CARTESIAN=cartesian, $
POLAR=polar, $
EXTEND=extend, $
ORDER=order, $
BASIS_ELLIPTICITY=basis_ellipticity, $
COEFFICIENTS=coefficients, $
NOHISTORY=nohistory, $
MAINTAIN=maintain
;$Id: shapelets_shear.pro, v2$
;
; Copyright © 2005 Richard Massey and Alexandre Refregier.
;
; This file is a part of the Shapelets analysis code.
; www.astro.caltech.edu/~rjm/shapelets/
;
; The Shapelets code is free software; you can redistribute it and/or
; modify it under the terms of the GNU General Public Licence as published
; by the Free Software Foundation; either version 2 of the Licence, or
; (at your option) any later version.
;
; The Shapelets code is distributed in the hope that it will be useful,
; but WITHOUT ANY WARRANTY; without even the implied warranty of
; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
; GNU General Public Licence for more details.
;
; You should have received a copy of the GNU General Public Licence
; along with the Shapelets code; if not, write to the Free Software
; Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
;
;+
; NAME:
; SHAPELETS_SHEAR
;
; CATEGORY:
; Shapelets.
;
; PURPOSE:
; Applies a (pure) shear to one object in a decomp structure, or to all
; objects in a shapecat catalogue.
;
; INPUTS:
; STRUCTURE - A shapelet decomp or shapecat structure.
; GAMMAS - Shear in [gamma_1,gamma_2] or complex(gamma1,gamma2) format.
;
; OPTIONAL INPUTS:
; EXTEND - Increses n_max before starting to shear, and set the new
; coefficients to zero. This gives the shear somewhere to go.
; ORDER - Include additional orders of shear operator by setting this
; to an integer higher than 1. Default is to apply only those
; terms of order gamma.
;
; KEYWORD PARAMETERS:
; BASIS_ELLI- Forces operation to be performed by changing the underlying
; ellipticity of the shapelet basis functions.
; COEFFICIEN- Forces operation to be performed via transformation of the
; shapelet coefficients.
; DEFAULT: Elliptical shapelet models are sheared by changing
; the ellipticity of the basis functions.
; Circular shapelets are sheared via coefficients.
; (be aware that two shears don't commute: see BJ02).
;
; OUTPUTS:
; STRUCTURE - The structure is returned, sheared.
;
; TO DO:
; Accept sexcats as an input. This would be an interesting exercise.
; Properly transform the errors on the shapelet coefficients.
;
; MODIFICATION HISTORY:
;;;;;;; Apr 09 - If order>1, include change of size but not change of flux RM.
; Jan 09 - If order>1, incorporate magnification from det(A) by RM.
; Jan 09 - If order>1, use BJ02 eta shear (rather than delta) by RM.
; May 07 - Handling of basis functions with theta_zero<>0 added by RM.
; Apr 07 - Handling of elliptical shapelet basis functions added by RM.
; Jul 05 - Sign error fixed in polar implementation by RM.
; Apr 05 - Shapecats accepted as input by RM.
; Apr 05 - POLAR keyword added by RM.
; Apr 05 - Higher order transformation added properly by RM.
; Feb 05 - IDL2 notation (square brackets, etc) compatability by RM.
; Aug 04 - Approximate higher order shear calculation added by RM.
; Apr 02 - Written by Richard Massey
;-
COMPILE_OPT idl2
; Maintain backwards compatibility
if not shapelets_structure_type(structure,message=message) then message,message
; Parse inputs
if size(gammas,/type) eq 6 then begin
gamma1=float(gammas[0])
gamma2=imaginary(gammas[0])
gammac=gammas[0]
endif else begin
gamma1=gammas[0]
if n_elements(gammas) ge 2 then gamma2=gammas[1] else gamma2=0
gammac=complex(gamma1,gamma2)
endelse
; Add operation to object's history record
if not keyword_set(nohistory) then shapelets_update_history, structure, $
"Sheared by {"+strmid(strtrim(gamma1,2),0,5)+","+strmid(strtrim(gamma2,2),0,5)+"}"
; Take quick route with objects that use elliptical shapelet basis functions
if tag_exist(structure,"basis_ellipticity") and keyword_set(basis_ellipticity) and not keyword_set(coefficients) then begin
if keyword_set(basis_ellipticity) or total(abs(structure.basis_ellipticity)) gt 0. then begin
print,'input basis ellipticity ',structure.basis_ellipticity
eta=abs(alog(1-abs(structure.basis_ellipticity))-alog(1+abs(structure.basis_ellipticity))) ; From BJ02 eq (2.8), with atanh(z)=[ln(1+z)-ln(1-z)]/2
print,'eta ',eta
theta=atan(imaginary(structure.basis_ellipticity),float(structure.basis_ellipticity))/2. ; but this includes the factor of 2
print,'theta',theta
delta1=tanh(eta)*[cos(2*theta),sin(2*theta)] ; From BJ02 eq (2.7)
message,/INFO,"Make sure elliptical shapelets are sheared by eta rather than delta!"
theta=atan(gamma2,gamma1)/2.
print,'theta',theta
delta2=tanh(atanh(abs(gammac))*2.)*[cos(2*theta),sin(2*theta)]
delta2sq=total(delta2^2)
delta3=[(delta1[0] + delta2[0] + $ ; BJ02 eq. (2.13)
(delta2[1]/delta2sq)*(1-sqrt(1-delta2sq))*(delta1[1]*delta2[0]-delta1[0]*delta2[1])) $
,(delta1[1] + delta2[1] + $
(delta2[0]/delta2sq)*(1-sqrt(1-delta2sq))*(delta1[0]*delta2[1]-delta1[1]*delta2[0])) ] $
/(1+delta1[0]*delta2[0]+delta1[1]*delta2[1])
print,'delta1',delta1
print,'delta2',delta2
print,'delta3',delta3
theta=atan(delta3[1],delta3[0])/2.
gamma3=tanh(atanh(sqrt(total(delta3^2)))/2.)*[cos(2*theta),sin(2*theta)]
print,'gamma3',gamma3
print,'theta',theta
;structure.beta=structure.beta*(1.+abs(structure.basis_ellipticity)^2)
structure.basis_ellipticity=complex(gamma3[0],gamma3[1])
print,'output basis ellipticity',structure.basis_ellipticity
print
;structure.beta=structure.beta/(1.+abs(structure.basis_ellipticity)^2)
return
endif
endif
; Parse more inputs
if not keyword_set(order) then order=1
if not keyword_set(cartesian) then cartesian=0B
if not keyword_set(polar) then polar=0B
polar_input=structure.polar
if tag_exist(structure,"theta_zero") then begin
if max(abs(structure.theta_zero)) gt 0 then begin
gammac=complex(gamma1*cos(2*structure.theta_zero/!radeg)+gamma2*sin(2*structure.theta_zero/!radeg),$
-gamma1*sin(2*structure.theta_zero/!radeg)+gamma2*cos(2*structure.theta_zero/!radeg))
gamma1=float(gammac[0])
gamma2=imaginary(gammac[0])
endif
endif
; Decide whether default method should be in Cartesian or polar shapelet space
if not ((cartesian+polar) mod 2) then begin
if structure.polar then begin
cartesian=0B
polar=1B
endif else begin
cartesian=1B
polar=0B
endelse
endif
; Perform shear to higher than first order
if fix(order) gt 1 then begin
shapelets_polar_convert,structure,CARTESIAN=cartesian,POLAR=polar,/SILENT
; Take into account the (second order) magnification by sqrt(1+gamma^2) during a shear - this is the determinant of the A matrix.
; However, do not change the flux (use BJ02 transformation)
shapelets_dilate,structure,1+abs(gammac)^2,/beta,/area,/SILENT;,/FLUX
; Use Bernstein & Jarvis (2002) eta shear rather than Miralda-Escude (1991) delta, where eta=arctanh(delta)
arctanh_coeffs=findgen(11>fix(order))+1 & arctanh_coeffs=[0,arctanh_coeffs mod 2/arctanh_coeffs]
shapelets_exponentiate_operations, "shapelets_shear", $
structure, poly(gammas,arctanh_coeffs), order, CARTESIAN=cartesian, POLAR=polar, EXTEND=extend
endif else begin
; Increase n_max to contain new coefficients that might be created
if keyword_set(extend) then shapelets_extend_nmax,structure,extend
; Perform shear to first order
if cartesian then begin
; Perform shear in Cartesian shapelet space
; Obtain initial shapelet coefficients
shapelets_polar_convert,structure,/CARTESIAN,/SILENT
if structure.type eq "shapecat" then begin
old_coeffs=transpose(structure.coeffs)
shapelets_make_nvec,structure.maxn_max,n1,n2,n_coeffs
endif else if structure.type eq "decomp" then begin
old_coeffs=structure.coeffs
shapelets_make_nvec,structure.n_max,n1,n2,n_coeffs
n1=structure.n1
n2=structure.n2
n_coeffs=structure.n_coeffs
endif else message,"Structure type not recognised!"
new_coeffs=old_coeffs
; Loop over all of the shapelet coefficients
for i=0,n_coeffs-1 do begin
; S_1
; Expand leftwards in Cartesian shapelet coefficient space
j=where(n1 eq n1[i]+2 and n2 eq n2[i],n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]-gamma1/2.*sqrt((n1[i]+1)*(n1[i]+2))*old_coeffs[j,*]
; Expand downwards in Cartesian shapelet coefficient space
j=where(n1 eq n1[i] and n2 eq n2[i]+2,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]+gamma1/2.*sqrt((n2[i]+1)*(n2[i]+2))*old_coeffs[j,*]
; Expand rightwards in Cartesian shapelet coefficient space
j=where(n1 eq n1[i]-2 and n2 eq n2[i],n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]+gamma1/2.*sqrt(n1[i]*(n1[i]-1))*old_coeffs[j,*]
; Expand upwards in Cartesian shapelet coefficient space
j=where(n1 eq n1[i] and n2 eq n2[i]-2,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]-gamma1/2.*sqrt(n2[i]*(n2[i]-1))*old_coeffs[j,*]
; S_2
; Expand diagonally down-left in Cartesian shapelet coefficient space
j=where(n1 eq n1[i]+1 and n2 eq n2[i]+1,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]-gamma2*sqrt((n1[i]+1)*(n2[i]+1))*old_coeffs[j,*]
; Expand diagonally up-right in Cartesian shapelet coefficient space
j=where(n1 eq n1[i]-1 and n2 eq n2[i]-1,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]+gamma2*sqrt(n1[i]*n2[i])*old_coeffs[j,*]
endfor
endif else begin
; Perform shear in polar shapelet space
; Obtain polar shapelet coefficients
shapelets_polar_convert,structure,/POLAR,/SILENT
if structure.type eq "shapecat" then begin
old_coeffs=transpose(structure.coeffs)
shapelets_make_nvec,structure.maxn_max,n,m,n_coeffs,/POLAR
endif else if structure.type eq "decomp" then begin
old_coeffs=structure.coeffs
;shapelets_make_nvec,structure.n_max,n,m,n_coeffs,/POLAR
n=structure.n
m=structure.m
n_coeffs=structure.n_coeffs
endif else message,"Structure type not recognised!"
new_coeffs=old_coeffs
; Loop over all of the shapelet coefficients
for i=0,n_coeffs-1 do begin
; Expand up-right in polar shapelet coefficient space
j=where(n eq n[i]-2 and m eq m[i]-2,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]+gammac/4.*sqrt((n[i]+m[i])*(n[i]+m[i]-2))*old_coeffs[j,*]
; Expand down-right in polar shapelet coefficient space
j=where(n eq n[i]+2 and m eq m[i]-2,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]-gammac/4.*sqrt((n[i]-m[i]+2)*(n[i]-m[i]+4))*old_coeffs[j,*]
; Expand up-left in polar shapelet coefficient space
j=where(n eq n[i]-2 and m eq m[i]+2,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]+conj(gammac)/4.*sqrt((n[i]-m[i])*(n[i]-m[i]-2))*old_coeffs[j,*]
; Expand down-left in polar shapelet coefficient space
j=where(n eq n[i]+2 and m eq m[i]+2,n_j)
if n_j gt 0 then new_coeffs[i,*]=new_coeffs[i,*]-conj(gammac)/4.*sqrt((n[i]+m[i]+2)*(n[i]+m[i]+4))*old_coeffs[j,*]
endfor
endelse
; Reinsert new coefficients into the old structure
if structure.type eq "shapecat" then begin
new_coeffs=transpose(new_coeffs)
if tag_exist(structure,"moments") then if structure.moments then $
message,"The moments in your shapecat need updating!",/info
endif
structure.coeffs=new_coeffs
endelse
; Recover input format
if keyword_set(maintain) then $
shapelets_polar_convert,structure,CARTESIAN=1-polar_input,POLAR=polar_input,/SILENT
end
Return to the shapelets web page or the code help menu.
| Last modified on 2nd Mar 2009 by Richard Massey. |