MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Hole/Disk function

  • To: mathgroup at smc.vnet.net
  • Subject: [mg76108] Re: [mg76024] Hole/Disk function
  • From: Daniel Lichtblau <danl at wolfram.com>
  • Date: Wed, 16 May 2007 05:35:33 -0400 (EDT)
  • References: <200705150845.EAA16651@smc.vnet.net>

Mathieu G wrote:
> Hello,
> Is there a better way to define a Disk/Hole function than that?:
> 
> DHole[x_, y_] := If[Sqrt[x^2 + y^2] <= HoleSize, 1, 0];
> 
> I would like to use it for a convolution using the definition as shown 
> in http://oldsite.vislab.usyd.edu.au/CP3/Four5/mod2-01.gif
> 
> Best regards,
> MG

I'd avoid If[] and instead use UnitStep, HeavisideTheta, Piecewise, or 
Boole to represent the cutoff. Also I suspect you want to normalize so 
that the integral over your aperture has unit weight. So you could do as 
below, for example.

holeSize = 150*^-9;
beamRadius = 381*^-4;

diskMask[x_, y_] :=
   Boole[x^2 + y^2 <= holeSize]/(Pi*holeSize^2)

squareMask[x_, y_] :=
   Boole[Abs[x]<=holeSize && Abs[y]<=holeSize]/(4*holeSize^2)

gaussian2D[x_, y_] := Exp[-1/2*((x/beamRadius)^2 + (y/beamRadius)^2)]

The symbolic form of convolution would be

convolve[f_, g_, x_, y_] :=  Integrate[
   g[xPrime, yPrime]*f[x - xPrime, y - yPrime], {xPrime, -Infinity,
    Infinity}, {yPrime, -Infinity, Infinity},
   GenerateConditions -> True]

You could instead evaluate numerically as follows.

convolveN[f_, g_, x_?NumericQ, y_?NumericQ] := NIntegrate[
   g[xPrime, yPrime]*f[x - xPrime, y - yPrime],
   {xPrime, -Infinity, Infinity}, {yPrime, -Infinity, Infinity}]

The integral for the square aperture is simple to obtain.

cdsquare[x_, y_] = convolve[squareMask, gaussian2D, x, y]

I'm not sure the circular case can be handled that way directly, though 
you can get numeric values as below.

cddiskN[x_?NumericQ, y_?NumericQ] =
   convolveN[squareMask, gaussian2D, x, y]

To obtain it symbolically you might try converting to polar coordinates 
and see if the convolution integral can be done. Alternatively, try to 
get it via Fourier transforms (or maybe combine these approaches).


Daniel Lichtblau
Wolfram Research






  • Prev by Date: Re: Solve & RotationMatrix
  • Next by Date: Re: weird behavior when plotting multiple functions in
  • Previous by thread: Re: Hole/Disk function
  • Next by thread: Re: Hole/Disk function