mperrin / poppy

Physical Optics Propagation in Python
BSD 3-Clause "New" or "Revised" License
177 stars 41 forks source link

Adding user-defined Fresnel wavefront error #261

Closed mb2448 closed 6 years ago

mb2448 commented 6 years ago

cf #1: Is it possible to add wavefront error to a Fresnel wavefront if it is not one of the standard wavefront error classes? The following code crashes:

import poppy
import numpy as np
import astropy.units as u

if __name__ == "__main__":
    wl = 500e-9
    #create a wavefront
    wf  = poppy.FresnelWavefront(0.01*u.m,
                                 wavelength=wl,
                                 npix=1024, 
                                 oversample = 2)
    #add 1 nm surface error through a "screen"
    surf_error = np.random.normal(loc = 0, scale = 1e-9, size = wf.shape)
    #convert to phase 
    error = np.exp(2*np.pi*1j/wl*surf_error)
    #next line crashes
    wf *= error

A scan through the API documentation did not suggest a solution.

mb2448 commented 6 years ago

Never mind, you can do this with the ArrayOpticalElement class, which I did not have with version 0.6. Here is an implementation below showing that it is doing the right thing

import poppy
import numpy as np
import astropy.units as u

if __name__ == "__main__":
    wl = 500e-9
    #create a wavefront
    wf  = poppy.FresnelWavefront(0.01*u.m,
                                 wavelength=wl,
                                 npix=1024,
                                 oversample = 2)
    #add 1 nm surface error through a "screen"
    surf_error = np.random.normal(loc = 0, scale = 1e-9, size = wf.shape)
    #convert to phase 
    error = np.exp(2*np.pi*1j/wl*surf_error)
    #Use the poppy version
    error_p= poppy.ArrayOpticalElement(opd = surf_error)
    #make sure the phasor is in accordance with expectation
    print np.allclose(error, error_p.get_phasor(wl*u.m))
    wf *= error_p
mperrin commented 6 years ago

Great, glad you figured this out, Michael.

The problem with the first version is that the *= invokes the Wavefront.__imul__ function, which knows how to "multiply" a wavefront by a scalar, an OpticalElement instance, or a coordinate transform, but not a plain ndarray. Wrapping the error into an ArrayOpticalElement fixes that, as you figured out.

The other way you could have done this in your first example is by directly accessing the wavefront data as an ndarray rather than via the Wavefront class:

wf.wavefront *= error

It would be simple enough to extend the Wavefront.__imul__ function to do the right thing if it's handed a complex ndarray with the right dimensions. Do you think that would be a useful enough thing to add in here?

Cleaning up and improving the documentation is on our to-do list for the next release this fall. I'll make sure we include an example along these lines!

mb2448 commented 6 years ago

Thanks Marshall. In my opinion it would be useful to allow for multiplication by complex matrices of the correct shape natively, without passing to an optical element class...but it seems straightforward enough to do it either way.