From 8c14f1d099ac990dd8dd47caf4574c859e0c1f09 Mon Sep 17 00:00:00 2001 From: A J Eldorado Riggs Date: Wed, 9 Jul 2025 17:51:20 -0500 Subject: [PATCH 1/3] AD-EFC: FIx cost function. Linting --- falco/ctrl.py | 75 +++++++++++++++++++++---------------------- falco/model/models.py | 2 +- 2 files changed, 37 insertions(+), 40 deletions(-) diff --git a/falco/ctrl.py b/falco/ctrl.py index a70fe36..ac1fc9e 100644 --- a/falco/ctrl.py +++ b/falco/ctrl.py @@ -44,17 +44,17 @@ def wrapper(mp, cvar, jacStruct): else: apply_spatial_weighting_to_Jacobian(mp, jacStruct) - + # with falco.util.TicToc('Using the Jacobian to make other matrices'): print('Using the Jacobian to make other matrices...', end='') - + # Compute matrices for linear control with regular EFC cvar.GstarG_wsum = np.zeros((cvar.NeleAll, cvar.NeleAll)) cvar.RealGstarEab_wsum = np.zeros((cvar.NeleAll, 1)) Eest = cvar.Eest.copy() - + for iMode in range(mp.jac.Nmode): - + Gstack = np.zeros((mp.Fend.corr.Npix, 1), dtype=complex) # Initialize a row to concatenate onto if any(mp.dm_ind == 1): Gstack = np.hstack((Gstack, np.squeeze(jacStruct.G1[:, :, iMode]))) @@ -65,10 +65,10 @@ def wrapper(mp, cvar, jacStruct): if any(mp.dm_ind == 9): Gstack = np.hstack((Gstack, np.squeeze(jacStruct.G9[:, :, iMode]))) Gstack = Gstack[:, 1:] # Remove the column used for initialization - + # Square matrix part stays the same if no re-linearization has occurrred. cvar.GstarG_wsum += mp.jac.weights[iMode]*np.real(np.conj(Gstack).T @ Gstack) - + if mp.jac.minimizeNI: modvar = falco.config.ModelVariables() modvar.whichSource = 'star' @@ -80,7 +80,7 @@ def wrapper(mp, cvar, jacStruct): np.argmax(np.abs(Eunocculted), axis=None), Eunocculted.shape) Epeak = Eunocculted[indPeak] Eest[:, iMode] = cvar.Eest[:, iMode] / Epeak - + # The G^*E part changes each iteration because the E-field changes. # Apply 2-D spatial weighting to E-field in dark hole pixels. iStar = mp.jac.star_inds[iMode] @@ -88,12 +88,12 @@ def wrapper(mp, cvar, jacStruct): # Apply the Jacobian weights and add to the total. cvar.RealGstarEab_wsum += mp.jac.weights[iMode]*np.real( np.conj(Gstack).T @ Eweighted.reshape(mp.Fend.corr.Npix, 1)) - + # Make the regularization matrix. (Define only diagonal here to save RAM.) cvar.EyeGstarGdiag = np.max(np.diag(cvar.GstarG_wsum))*np.ones(cvar.NeleAll) cvar.EyeNorm = np.max(np.diag(cvar.GstarG_wsum)) print('done.') - + # Call the Controller Function print('Control beginning ...') # Established, conventional controllers @@ -448,7 +448,7 @@ def _planned_ad_efc(mp, cvar): # Step 1: Empirically find the "optimal" regularization value # (if told to for this iteration). - + if runNewGridSearch: # Temporarily store computed DM commands so that the best one does # not have to be re-computed @@ -460,14 +460,14 @@ def _planned_ad_efc(mp, cvar): dDM8V_store = np.zeros((mp.dm8.NactTotal, Nvals)) if any(mp.dm_ind == 9): dDM9V_store = np.zeros((mp.dm9.NactTotal, Nvals)) - + ImCube = np.zeros((Nvals, mp.Fend.Neta, mp.Fend.Nxi)) - + for ni in range(Nvals): - + [InormVec[ni], dDM_temp] = _ad_efc(ni, vals_list, mp, cvar) ImCube[ni, :, :] = dDM_temp.Itotal - + # delta voltage commands if any(mp.dm_ind == 1): dDM1V_store[:, :, ni] = dDM_temp.dDM1V @@ -477,21 +477,21 @@ def _planned_ad_efc(mp, cvar): dDM8V_store[:, ni] = dDM_temp.dDM8V if any(mp.dm_ind == 9): dDM9V_store[:, ni] = dDM_temp.dDM9V - + # Print out results to the command line print('Scaling factor:\t\t', end='') for ni in range(Nvals): print('%.2f\t\t' % (vals_list[ni][1]), end='') - + print('\nlog10reg: \t\t', end='') for ni in range(Nvals): print('%.1f\t\t' % (vals_list[ni][0]), end='') - + print('\nInorm: \t\t', end='') for ni in range(Nvals): print('%.2e\t' % (InormVec[ni]), end='') print('\n', end='') - + # Find the best scaling factor and Lagrange multiplier pair based on # the best contrast. # [cvar.cMin,indBest] = np.min(InormVec) @@ -500,7 +500,7 @@ def _planned_ad_efc(mp, cvar): cvar.Im = np.squeeze(ImCube[indBest, :, :]) cvar.latestBestlog10reg = vals_list[indBest][0] cvar.latestBestDMfac = vals_list[indBest][1] - + if mp.ctrl.flagUseModel: print(('Model-based grid search expects log10reg, = %.1f,\t ' + 'dmfac = %.2f,\t %4.2e normalized intensity.') % @@ -509,7 +509,7 @@ def _planned_ad_efc(mp, cvar): print(('Empirical grid search finds log10reg, = %.1f,\t dmfac' + ' = %.2f,\t %4.2e normalized intensity.') % (cvar.latestBestlog10reg, cvar.latestBestDMfac, cvar.cMin)) - + # Skip steps 2 and 3 if the schedule for this iteration is just to use the # "optimal" regularization AND if grid search was performed this iteration. if runNewGridSearch and useBestLog10Reg and realLog10RegIsZero: @@ -523,7 +523,7 @@ def _planned_ad_efc(mp, cvar): dDM.dDM8V = np.squeeze(dDM8V_store[:, indBest]) if any(mp.dm_ind == 9): dDM.dDM9V = np.squeeze(dDM9V_store[:, indBest]) - + log10regSchedOut = cvar.latestBestlog10reg else: # Step 2: For this iteration in the schedule, replace the imaginary @@ -875,8 +875,8 @@ def _efc(ni, vals_list, mp, cvar): def _ad_efc(ni, vals_list, mp, cvar): """ - Use algorithmic differentiation to suppress the E-field. - + Use algorithmic differentiation to suppress the E-field. + Called by a wrapper controller function. Parameters @@ -967,8 +967,6 @@ def _ad_efc(ni, vals_list, mp, cvar): EFend = falco.model.compact(mp, modvar, isNorm=True, isEvalMode=False, useFPM=True, forRevGradModel=False) EFendPrev.append(EFend) - - t0 = time.time() @@ -977,7 +975,7 @@ def _ad_efc(ni, vals_list, mp, cvar): method='L-BFGS-B', jac=True, bounds=bounds, tol=None, callback=None, options={'disp': None, - 'ftol': 1e-12, + 'ftol': 1e-99, 'gtol': 1e-10, 'maxls': 20, 'maxiter': mp.ctrl.ad.maxiter, @@ -989,19 +987,18 @@ def _ad_efc(ni, vals_list, mp, cvar): duVec = u_sol.x print(u_sol.success) print(u_sol.nit) - + ''' def optim_cost_func(Vdm): - #cost function for L-BFGS + #cost function for L-BFGS return falco.model.compact_reverse_gradient(Vdm, mp, cvar.Eest, EFendPrev, log10reg) - + opt = F77LBFGSB(fg=optim_cost_func, x0=dm0) nIter = 0; for _ in range(mp.ctrl.ad.maxiter): nIter += 1 xk, fk, gk = opt.step() - - + print("Niter LBFGS: ", nIter) duVec = xk ''' @@ -1102,18 +1099,18 @@ def init(mp, cvar): cvar.DM9Vnom = mp.dm9.V # Make the vector of total DM commands from before - u1 = mp.dm1.V.reshape(mp.dm1.NactTotal)[mp.dm1.act_ele] if(any(mp.dm_ind == 1)) else np.array([]) - u2 = mp.dm2.V.reshape(mp.dm2.NactTotal)[mp.dm2.act_ele] if(any(mp.dm_ind == 2)) else np.array([]) - u8 = mp.dm9.V.reshape(mp.dm8.NactTotal)[mp.dm8.act_ele] if(any(mp.dm_ind == 8)) else np.array([]) - u9 = mp.dm9.V.reshape(mp.dm9.NactTotal)[mp.dm9.act_ele] if(any(mp.dm_ind == 9)) else np.array([]) + u1 = mp.dm1.V.reshape(mp.dm1.NactTotal)[mp.dm1.act_ele] if any(mp.dm_ind == 1) else np.array([]) + u2 = mp.dm2.V.reshape(mp.dm2.NactTotal)[mp.dm2.act_ele] if any(mp.dm_ind == 2) else np.array([]) + u8 = mp.dm9.V.reshape(mp.dm8.NactTotal)[mp.dm8.act_ele] if any(mp.dm_ind == 8) else np.array([]) + u9 = mp.dm9.V.reshape(mp.dm9.NactTotal)[mp.dm9.act_ele] if any(mp.dm_ind == 9) else np.array([]) cvar.uVec = np.concatenate((u1, u2, u8, u9)) cvar.NeleAll = cvar.uVec.size # Get the indices of each DM's command within the full command - u1dummy = 1*np.ones(mp.dm1.Nele, dtype=int) if(any(mp.dm_ind == 1)) else np.array([]) - u2dummy = 2*np.ones(mp.dm2.Nele, dtype=int) if(any(mp.dm_ind == 2)) else np.array([]) - u8dummy = 8*np.ones(mp.dm8.Nele, dtype=int) if(any(mp.dm_ind == 8)) else np.array([]) - u9dummy = 9*np.ones(mp.dm9.Nele, dtype=int) if(any(mp.dm_ind == 9)) else np.array([]) + u1dummy = 1*np.ones(mp.dm1.Nele, dtype=int) if any(mp.dm_ind == 1) else np.array([]) + u2dummy = 2*np.ones(mp.dm2.Nele, dtype=int) if any(mp.dm_ind == 2) else np.array([]) + u8dummy = 8*np.ones(mp.dm8.Nele, dtype=int) if any(mp.dm_ind == 8) else np.array([]) + u9dummy = 9*np.ones(mp.dm9.Nele, dtype=int) if any(mp.dm_ind == 9) else np.array([]) cvar.uLegend = np.concatenate((u1dummy, u2dummy, u8dummy, u9dummy)) mp.ctrl.uLegend = cvar.uLegend diff --git a/falco/model/models.py b/falco/model/models.py index d21f9a6..d66e20a 100644 --- a/falco/model/models.py +++ b/falco/model/models.py @@ -537,7 +537,7 @@ def compact_reverse_gradient(command_vec, mp, EestAll, EFendPrev, log10reg): DH = EdhNew[mp.Fend.corr.maskBool] int_in_dh = np.sum(np.abs(DH)**2) # print(f'int_in_dh = {int_in_dh}') - total_cost += mp.jac.weights[imode] * int_in_dh / normFacAD + total_cost += mp.jac.weights[imode] * int_in_dh normFacADweightedSum += mp.jac.weights[imode] / normFacAD # Gradient From a65527aea66928a19a1d79ffc983f310c0279c5a Mon Sep 17 00:00:00 2001 From: A J Eldorado Riggs Date: Wed, 9 Jul 2025 17:57:54 -0500 Subject: [PATCH 2/3] Make PTP propagation return input if zero propagation distance --- falco/prop.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/falco/prop.py b/falco/prop.py index ab55c15..1b91514 100644 --- a/falco/prop.py +++ b/falco/prop.py @@ -83,6 +83,10 @@ def ptp(E_in, full_width, wavelength, dz): check.real_positive_scalar(wavelength, 'wavelength', TypeError) check.real_scalar(dz, 'dz', TypeError) + # Return unchanged if no propagation + if dz == 0: + return E_in + M, N = E_in.shape dx = full_width / N N_critical = int(np.floor(wavelength * np.abs(dz) / (dx ** 2))) From e6538339bb44a6ed45fab10f6e103e80d6ced1bd Mon Sep 17 00:00:00 2001 From: A J Eldorado Riggs Date: Wed, 9 Jul 2025 18:20:22 -0500 Subject: [PATCH 3/3] Fix bug where single DM usage was erroring And add a functional test to prove that single DM EFC works. --- falco/model/jacobians.py | 5 +- .../functional/config_wfsc_lc_quick_one_dm.py | 297 ++++++++++++++++++ tests/functional/test_wfsc_lc_one_dm.py | 63 ++++ 3 files changed, 363 insertions(+), 2 deletions(-) create mode 100644 tests/functional/config_wfsc_lc_quick_one_dm.py create mode 100644 tests/functional/test_wfsc_lc_one_dm.py diff --git a/falco/model/jacobians.py b/falco/model/jacobians.py index ffb09f9..5d5f115 100644 --- a/falco/model/jacobians.py +++ b/falco/model/jacobians.py @@ -186,8 +186,9 @@ def precomp(mp): Edm1 *= Edm1WFE*DM1stop*np.exp(surfIntoPhase*2*np.pi*1j*DM1surf/wvl) # Two array sizes (at same resolution) of influence functions for MFT and angular spectrum - NboxPad2AS = int(mp.dm2.compact.NboxAS) - mp.dm2.compact.xy_box_lowerLeft_AS = mp.dm2.compact.xy_box_lowerLeft - (NboxPad2AS-mp.dm2.compact.Nbox)/2 # Account for the padding of the influence function boxes + if 2 in mp.dm_ind: + NboxPad2AS = int(mp.dm2.compact.NboxAS) + mp.dm2.compact.xy_box_lowerLeft_AS = mp.dm2.compact.xy_box_lowerLeft - (NboxPad2AS-mp.dm2.compact.Nbox)/2 # Account for the padding of the influence function boxes apodReimaged = pad_crop(apodReimaged, mp.dm2.compact.NdmPad) DM2stopPad = pad_crop(DM2stop, mp.dm2.compact.NdmPad) diff --git a/tests/functional/config_wfsc_lc_quick_one_dm.py b/tests/functional/config_wfsc_lc_quick_one_dm.py new file mode 100644 index 0000000..fd0083a --- /dev/null +++ b/tests/functional/config_wfsc_lc_quick_one_dm.py @@ -0,0 +1,297 @@ +"""Script to define configuration values for FALCO.""" +import numpy as np + +import falco + +mp = falco.config.ModelParameters() +mp.path = falco.config.Object() + +# Record Keeping +mp.SeriesNum = 1 +mp.TrialNum = 1 + +# Special Computational Settings +mp.flagParallel = False +mp.flagPlot = False + +# General +mp.centering = 'pixel' + +# Method of computing core throughput: +# - 'HMI' for energy within half-max isophote divided by energy at telescope pupil +# - 'EE' for encircled energy within a radius (mp.thput_radius) divided by energy at telescope pupil +mp.thput_metric = 'HMI' +mp.thput_radius = 0.7 # photometric aperture radius [lambda_c/D]. Used ONLY for 'EE' method. +mp.thput_eval_x = 7 # x location [lambda_c/D] in dark hole at which to evaluate throughput +mp.thput_eval_y = 0 # y location [lambda_c/D] in dark hole at which to evaluate throughput + +# Where to shift the source to compute the intensity normalization value. +mp.source_x_offset_norm = 7 # x location [lambda_c/D] in dark hole at which to compute intensity normalization +mp.source_y_offset_norm = 0 # y location [lambda_c/D] in dark hole at which to compute intensity normalization + +# Bandwidth and Wavelength Specs +mp.lambda0 = 575e-9 # Central wavelength of the whole spectral bandpass [meters] +mp.fracBW = 0.01 # fractional bandwidth of the whole bandpass (Delta lambda / lambda0) +mp.Nsbp = 1 # Number of sub-bandpasses to divide the whole bandpass into for estimation and control +mp.Nwpsbp = 1 # Number of wavelengths to used to approximate an image in each sub-bandpass + + +# %% Wavefront Estimation + +# Estimator Options: +# - 'perfect' for exact numerical answer from full model +# - 'pwp-bp' for pairwise probing with batch process estimation +mp.estimator = 'perfect' + +# Pairwise probing: +mp.est = falco.config.Object() +mp.est.probe = falco.config.Probe() +mp.est.probe.Npairs = 3 # Number of pair-wise probe PAIRS to use. +mp.est.probe.whichDM = 1 # Which DM # to use for probing. 1 or 2. Default is 1 +mp.est.probe.radius = 12 # Max x/y extent of probed region [lambda/D]. +mp.est.probe.xOffset = 0 # offset of probe center in x [lambda/D]. Use to avoid central obscurations. +mp.est.probe.yOffset = 14 # offset of probe center in y [lambda/D]. Use to avoid central obscurations. +mp.est.probe.axis = 'alternate' # which axis to have the phase discontinuity along [x or y or xy/alt/alternate] +mp.est.probe.gainFudge = 1 # empirical fudge factor to make average probe amplitude match desired value. + + +# %% Wavefront Control: General + +mp.ctrl = falco.config.Object() +mp.ctrl.flagUseModel = False # Whether to perform a model-based (vs empirical) grid search for the controller + +# Threshold for culling weak actuators from the Jacobian: +mp.logGmin = -6 # 10^(mp.logGmin) used on the intensity of DM1 and DM2 Jacobians to weed out the weakest actuators + +# Zernikes to suppress with controller +mp.jac = falco.config.Object() +mp.jac.zerns = [1, ] # Which Noll Zernike modes to include in Jacobian. Always include the value "1" for the on-axis piston mode. +mp.jac.Zcoef = 1e-9*np.ones(np.size(mp.jac.zerns)) # meters RMS of Zernike aberrations. (piston value is reset to 1 later) + +# Zernikes to compute sensitivities for +mp.eval = falco.config.Object() +mp.eval.indsZnoll = [2, 3] # Noll indices of Zernikes to compute values for [1-D ndarray] + +# Annuli to compute 1nm RMS Zernike sensitivities over. Columns are [inner radius, outer radius]. One row per annulus. +mp.eval.Rsens = np.array([[3., 4.], [4., 8.]]) # [2-D ndarray] + +# Grid- or Line-Search Settings +mp.ctrl.log10regVec = np.arange(-6, -2+0.5, 1) # log10 of the regularization exponents (often called Beta values) +mp.ctrl.dmfacVec = np.array([1.]) # Proportional gain term applied to the total DM delta command. Usually in range [0.5,1]. [1-D ndarray] + +# Spatial pixel weighting +mp.WspatialDef = [] # [3, 4.5, 3]; #--spatial control Jacobian weighting by annulus: [Inner radius, outer radius, intensity weight; (as many rows as desired)] [ndarray] + +# DM weighting +mp.dm1.weight = 1. +mp.dm2.weight = 1. + + +# %% Wavefront Control: Controller Specific +# Controller options: +# - 'gridsearchEFC' for EFC as an empirical grid search over tuning parameters +# - 'plannedEFC' for EFC with an automated regularization schedule +mp.controller = 'gridsearchEFC' + +# # GRID SEARCH EFC DEFAULTS +# WFSC Iterations and Control Matrix Relinearization +mp.Nitr = 3 # Number of estimation+control iterations to perform +mp.relinItrVec = np.arange(mp.Nitr) # Which correction iterations at which to re-compute the control Jacobian [1-D ndarray] +mp.dm_ind = [1, ] # Which DMs to use [1-D array_like] + +# # PLANNED SEARCH EFC DEFAULTS +# mp.dm_ind = [1 2 ]; # vector of DMs used in controller at ANY time (not necessarily all at once or all the time). +# mp.ctrl.dmfacVec = 1; +# #--CONTROL SCHEDULE. Columns of mp.ctrl.sched_mat are: +# # Column 1: # of iterations, +# # Column 2: log10(regularization), +# # Column 3: which DMs to use (12, 128, 129, or 1289) for control +# # Column 4: flag (0 = False, 1 = True), whether to re-linearize +# # at that iteration. +# # Column 5: flag (0 = False, 1 = True), whether to perform an +# # EFC parameter grid search to find the set giving the best +# # contrast . +# # The imaginary part of the log10(regularization) in column 2 is +# # replaced for that iteration with the optimal log10(regularization) +# # A row starting with [0, 0, 0, 1...] is for relinearizing only at that time +# +# mp.ctrl.sched_mat = [... +# repmat([1,1j,12,1,1],[4,1]);... +# repmat([1,1j-1,12,1,1],[25,1]);... +# repmat([1,1j,12,1,1],[1,1]);... +# ]; +# [mp.Nitr, mp.relinItrVec, mp.gridSearchItrVec, mp.ctrl.log10regSchedIn, mp.dm_ind_sched] = falco_ctrl_EFC_schedule_generator(mp.ctrl.sched_mat); + + +# %% Deformable Mirrors: Influence Functions + +# Influence Function Options: +# - falco.INFLUENCE_XINETICS uses the file 'influence_dm5v2.fits' for one type of Xinetics DM +# - INFLUENCE_BMC_2K uses the file 'influence_BMC_2kDM_400micron_res10.fits' for BMC 2k DM +# - INFLUENCE_BMC_KILO uses the file 'influence_BMC_kiloDM_300micron_res10_spline.fits' for BMC kiloDM +# mp.dm1.inf_fn = falco.INFLUENCE_BMC_2K +# mp.dm2.inf_fn = falco.INFLUENCE_BMC_2K +mp.dm1.inf_fn = falco.INFLUENCE_XINETICS +mp.dm2.inf_fn = falco.INFLUENCE_XINETICS + +mp.dm1.dm_spacing = 0.9906e-3 # actuator pitch +mp.dm2.dm_spacing = 0.9906e-3 # actuator pitch + +mp.dm1.inf_sign = '+' +mp.dm2.inf_sign = '+' + +# %% Deformable Mirrors: Optical Layout Parameters + +# DM1 parameters +mp.dm1.Nact = 32 # of actuators across DM array +mp.dm1.VtoH = 1e-9*np.ones((32, 32)) # gains of all actuators [nm/V of free stroke] +mp.dm1.xtilt = 0 # for foreshortening. angle of rotation about x-axis [degrees] +mp.dm1.ytilt = 5.83 # for foreshortening. angle of rotation about y-axis [degrees] +mp.dm1.zrot = 0 # clocking of DM surface [degrees] +mp.dm1.xc = (32/2 - 1/2) # x-center location of DM surface [actuator widths] +mp.dm1.yc = (32/2 - 1/2) # y-center location of DM surface [actuator widths] +mp.dm1.edgeBuffer = 1 # max radius (in actuator spacings) outside of beam on DM surface to compute influence functions for. [actuator widths] + +# DM2 parameters +mp.dm2.Nact = 32 # # of actuators across DM array +mp.dm2.VtoH = 1e-9*np.ones((32, 32)) # gains of all actuators [nm/V of free stroke] +mp.dm2.xtilt = 0 # for foreshortening. angle of rotation about x-axis [degrees] +mp.dm2.ytilt = 5.55 # for foreshortening. angle of rotation about y-axis [degrees] +mp.dm2.zrot = 0 # clocking of DM surface [degrees] +mp.dm2.xc = (32/2 - 1/2) # x-center location of DM surface [actuator widths] +mp.dm2.yc = (32/2 - 1/2) # y-center location of DM surface [actuator widths] +mp.dm2.edgeBuffer = 1 # max radius (in actuator spacings) outside of beam on DM surface to compute influence functions for. [actuator widths] + +# Aperture stops at DMs +mp.flagDM1stop = False # Whether to apply an iris or not +mp.dm1.Dstop = 50e-3 # Diameter of iris [meters] +mp.flagDM2stop = False # Whether to apply an iris or not +mp.dm2.Dstop = 50e-3 # Diameter of iris [meters] + +# DM separations +mp.d_P2_dm1 = 0 # distance (along +z axis) from P2 pupil to DM1 [meters] +mp.d_dm1_dm2 = 0 # distance between DM1 and DM2 [meters] + + +# %% Optical Layout: All models + +# Key Optical Layout Choices +mp.flagSim = True # Simulation or not +mp.layout = 'Fourier' # Which optical layout to use +mp.coro = 'LC' + +mp.Fend = falco.config.Object() + +# Final Focal Plane Properties +mp.Fend.res = 2.5 # Sampling [ pixels per lambda0/D] +mp.Fend.FOV = 11. # half-width of the field of view in both dimensions [lambda0/D] + +# Correction and scoring region definition +mp.Fend.corr = falco.config.Object() +mp.Fend.corr.Rin = 1.5 # inner radius of dark hole correction region [lambda0/D] +mp.Fend.corr.Rout = 9.7 # outer radius of dark hole correction region [lambda0/D] +mp.Fend.corr.ang = 180 # angular opening of dark hole correction region [degrees] +# +mp.Fend.score = falco.config.Object() +mp.Fend.score.Rin = 3.0 # inner radius of dark hole scoring region [lambda0/D] +mp.Fend.score.Rout = 9 # outer radius of dark hole scoring region [lambda0/D] +mp.Fend.score.ang = 180 # angular opening of dark hole scoring region [degrees] + +mp.Fend.sides = 'right' # Which side(s) for correction: 'left', 'right', 'top', 'up', 'bottom', 'down', 'lr', 'rl', 'leftright', 'rightleft', 'tb', 'bt', 'ud', 'du', 'topbottom', 'bottomtop', 'updown', 'downup' + +# %% Optical Layout: Compact Model (and Jacobian Model) + +# Focal Lengths +mp.fl = 1. # [meters] Focal length value used for all FTs in the compact model. Don't need different values since this is a Fourier model. + +# Pupil Plane Diameters +mp.P2.D = 30e-3 # [meters] +mp.P3.D = mp.P2.D # [meters] +mp.P4.D = mp.P2.D # [meters] + +# Pupil Plane Resolutions +mp.P1.compact.Nbeam = 100 +# mp.P2.compact.Nbeam = mp.P1.compact.Nbeam +# mp.P3.compact.Nbeam = mp.P1.compact.Nbeam +mp.P4.compact.Nbeam = mp.P1.compact.Nbeam # P4 size must be the same as P1 for Vortex. + +# Number of re-imaging relays between pupil planesin compact model. Needed +# to keep track of 180-degree rotations compared to the full model, which +# in general can have probably has extra collimated beams compared to the +# compact model. +mp.flagRotation = True # Whether to rotate 180 degrees between conjugate planes in the compact model +mp.Nrelay1to2 = 1 +mp.Nrelay2to3 = 1 +mp.Nrelay3to4 = 1 +mp.NrelayFend = 0 # How many times to rotate the final image by 180 degrees + +mp.F3.compact.res = 3 # sampling of FPM for compact model [pixels per lambda0/D] + +# %% Optical Layout: Full Model + +# # Focal Lengths +# mp.fl = 1 + +# Pupil Plane Resolutions +mp.P1.full.Nbeam = mp.P1.compact.Nbeam +# mp.P2.full.Nbeam = 250 +# mp.P3.full.Nbeam = 250 +mp.P4.full.Nbeam = mp.P1.full.Nbeam # P4 size must be the same as P1 for Vortex. + +# Mask Definitions +mp.full = falco.config.Object() +mp.compact = falco.config.Object() + +mp.F3.full.res = 4 # sampling of FPM for full model [pixels per lambda0/D] + + +# %% Entrance Pupil (P1) Definition and Generation + +mp.whichPupil = 'ROMAN' # Used only for run label +mp.P1.IDnorm = 0.303 # ID of the central obscuration [diameter]. Used only for computing the RMS DM surface from the ID to the OD of the pupil. OD is assumed to be 1. +mp.P1.ODnorm = 1.00 # Outer diameter of the telescope [diameter] +mp.P1.D = 2.3631 # circumscribed telescope diameter [meters]. Used only for converting milliarcseconds to lambda0/D or vice-versa. + +mp.P1.full.mask = falco.mask.falco_gen_pupil_Roman_CGI_20200513(mp.P1.full.Nbeam, mp.centering) +mp.P1.compact.mask = falco.mask.falco_gen_pupil_Roman_CGI_20200513(mp.P1.compact.Nbeam, mp.centering) + + +# %% "Apodizer" (P3) Definition and Generation + +mp.flagApod = False # Whether to use an apodizer or not + + +# %% Lyot stop (P4) Definition and Generation + +changes = {} +changes['flagLyot'] = True +changes['ID'] = 0.50 +changes['OD'] = 0.80 +changes['wStrut'] = 3.6/100 # nominal pupil's value is 76mm = 3.216% +changes['flagRot180'] = True + +mp.P4.full.mask = falco.mask.falco_gen_pupil_Roman_CGI_20200513(mp.P4.full.Nbeam, mp.centering, changes=changes) +mp.P4.compact.mask = falco.mask.falco_gen_pupil_Roman_CGI_20200513(mp.P4.compact.Nbeam, mp.centering, changes=changes) +mp.P4.compact.maskAtP1res = falco.mask.falco_gen_pupil_Roman_CGI_20200513(mp.P1.compact.Nbeam, mp.centering, changes=changes) + + +# %% FPM (F3) Definition and Generation + +mp.F3.Rin = 2.7 # maximum radius of inner part of the focal plane mask [lambda0/D] +mp.F3.RinA = 2.7 # inner hard-edge radius of the focal plane mask [lambda0/D]. Needs to be <= mp.F3.Rin +mp.F3.Rout = np.inf # radius of outer opaque edge of FPM [lambda0/D] +mp.F3.ang = 180 # on each side, opening angle [degrees] +mp.FPMampFac = 0 # amplitude transmission of the FPM + +# Both models +fpmDict = {} +fpmDict['rhoInner'] = mp.F3.Rin # radius of inner FPM amplitude spot (in lambda_c/D) +fpmDict['rhoOuter'] = mp.F3.Rout # radius of outer opaque FPM ring (in lambda_c/D) +fpmDict['centering'] = mp.centering +fpmDict['FPMampFac'] = mp.FPMampFac # amplitude transmission of inner FPM spot +# Full model +fpmDict['pixresFPM'] = mp.F3.full.res +mp.F3.full.mask = falco.mask.gen_annular_fpm(fpmDict) +# Compact model +fpmDict['pixresFPM'] = mp.F3.compact.res +mp.F3.compact.mask = falco.mask.gen_annular_fpm(fpmDict) diff --git a/tests/functional/test_wfsc_lc_one_dm.py b/tests/functional/test_wfsc_lc_one_dm.py new file mode 100644 index 0000000..9a8b3d9 --- /dev/null +++ b/tests/functional/test_wfsc_lc_one_dm.py @@ -0,0 +1,63 @@ +"""FALCO regression test of WFSC with a Lyot coronagraph.""" +from copy import deepcopy +import numpy as np +import os +from math import isclose + +import falco + +import config_wfsc_lc_quick_one_dm as CONFIG + + +def test_wfsc_lc(): + + mp = deepcopy(CONFIG.mp) + mp.flagPlot = False + LOCAL_PATH = os.path.dirname(os.path.abspath(__file__)) + mp.path.falco = os.path.dirname(os.path.dirname(LOCAL_PATH)) + + # Generate the label associated with this trial + mp.runLabel = 'testing_wfsc_lc' + + # Perform the Wavefront Sensing and Control + out = falco.setup.flesh_out_workspace(mp) + falco.wfsc.loop(mp, out) + + # print(out.IrawScoreHist[-1]) + + # print(out.IrawCorrHist[-1]) + # print(out.IestScoreHist[-1]) + # print(out.dm1.Spv[-1]) + # print(out.thput[-1]) + # print(out.log10regHist) + + # Tests: + print(f'out.IrawCorrHist[-1] = {out.IrawCorrHist[-1]}') + print(f'out.IestScoreHist[-1] = {out.IestScoreHist[-1]}') + print(f'out.dm1[-1] = {out.dm1.Spv[-1]}') + print(f'out.thput[-1] = {out.thput[-1]}') + print(f'out.log10regHist = {out.log10regHist}') + # out.IrawCorrHist[-1] = 5.049915478269254e-07 + # out.IestScoreHist[-1] = 8.070818118695458e-07 + # out.dm1[-1] = 7.684291769284934e-08 + # out.thput[-1] = 0.0768364858394983 + # out.log10regHist = [-2. -3. -3.] + + + Iend = out.IrawCorrHist[-1] + assert isclose(Iend, 5.049915478269254e-07, abs_tol=1e-8) + + Iest = out.IestScoreHist[-1] + assert isclose(Iest, 8.070818118695458e-07, abs_tol=1e-8) + + dm1pv = out.dm1.Spv[-1] + assert isclose(dm1pv, 7.684291769284934e-08, abs_tol=1e-9) + + thput = out.thput[-1] + assert isclose(thput, 0.0768364858394983, abs_tol=1e-3) + + assert np.allclose(out.log10regHist, np.array([-2, -3, -3]), rtol=1e-2) + + +if __name__ == '__main__': + test_wfsc_lc()