rconan / CEO

Cuda-Engined Adaptive Optics
zlib License
24 stars 18 forks source link

Glitch in the GMT pupil #55

Closed mvandam closed 3 years ago

mvandam commented 3 years ago

Below is an example of where the GMT pupil has glitches around the edges. This problem goes away for an odd number of pixels.

import ceo
import sys
import datetime
import os.path
import numpy as np
import matplotlib.pyplot as plt
plt.ion()

gmt = ceo.GMT_MX()       

D = 25.5
gmt.M2_baffle = 3.5    # Circular diameter pupil mask obscuration
gmt.project_truss_onaxis = True

nPx = 512 # number of pixels across D 
gs = ceo.Source("R+I",zenith=0.,azimuth=0.,rays_box_size=D, rays_box_sampling=nPx, rays_origin=[0.0,0.0,25])

gs>>(gmt,)

~gmt
+gs
amplitude = gs.wavefront.amplitude.host().ravel()
phase_ref = gs.wavefront.phase.host().ravel()
validpts = np.where(amplitude)[0]
nphasevals = len(validpts)

~gmt
state = gmt.state
state['M1']['Rxyz'][3,1] = 1e-6
gmt^=state
+gs

temp = np.full(nPx**2, np.nan)
temp[validpts] = (gs.wavefront.phase.host().ravel() - phase_ref)[validpts]
plt.figure(1)
plt.clf()
wf = np.reshape(temp,(nPx,nPx))
plt.imshow(wf,interpolation='None', origin='lower')
plt.colorbar()
plt.show()

print(wf[80,337])
print(wf[80,338])
print(wf[80,339]) # this pixel is wrong!
print()
print(wf[84,337])
print(wf[84,338])
print(wf[84,339]) # this pixel is wrong!
rconan commented 3 years ago

Yes this is a vignetting issue, the vignetting by the segments is recomputed for each propagation (+gs) and the vignetting is a function of the position of the segments (state) so you need to take into account the amplitude before and after applying the state when you differentiate the wavefront, like so:

import ceo
import sys
import datetime
import os.path
import numpy as np
import matplotlib.pyplot as plt
plt.ion()

gmt = ceo.GMT_MX()

D = 25.5
gmt.M2_baffle = 3.5 # Circular diameter pupil mask obscuration
gmt.project_truss_onaxis = True

nPx = 512 # number of pixels across D
gs = ceo.Source("R+I",zenith=0.,azimuth=0.,rays_box_size=D, rays_box_sampling=nPx, rays_origin=[0.0,0.0,25])

gs>>(gmt,)

~gmt
+gs
amplitude = gs.wavefront.amplitude.host().ravel()
phase_ref = gs.wavefront.phase.host().ravel()

~gmt
state = gmt.state
state['M1']['Rxyz'][3,1] = 1e-6
gmt^=state
+gs

# This is the only modification needed to the code 
validpts = np.where(amplitude*gs.wavefront.amplitude.host().ravel())[0]
nphasevals = len(validpts)

temp = np.full(nPx**2, np.nan)
temp[validpts] = (gs.wavefront.phase.host().ravel() - phase_ref)[validpts]
plt.figure(1)
plt.clf()
wf = np.reshape(temp,(nPx,nPx))
plt.imshow(wf,interpolation='None', origin='lower')
plt.colorbar()
plt.show()

print(wf[80,337])
print(wf[80,338])
print(wf[80,339]) # this pixel is [wrong] vignetted!
print()
print(wf[84,337])
print(wf[84,338])
print(wf[84,339]) # this pixel is [wrong] vignetted!