基于MATLAB的GS算法

%
% GS_2D.m
%
% Calculation of a 2D phase or amplitude hologram by using the Gerchberg-Saxton
% algorithm.
%
% The program starts with a plan wave （a phase hologram with all pixels set
% to the same value） （i.e an electric field） and by iteration the operation
% of this hologram is increased.
%
% The program use the FFT-file fft2c.m and the inverse IFFT-file ifft2c.m
%
% All lengths are given in micrometers.
%
%
% Orginal file was created during the spring of 2001 by David Engstr鰉 Chalmers
% 
% Modified in November 2002 by Emil H鋖lstig FOI
%

close all; % Close all figures
clear all; % Clear workspace

% Constants in use:


N=256;	    % The number of sample points in one direction （totally N*N sampling points N
            % should be set to 2^integer to increase the speed of the FFT-transform）

nr_of_iterations=10;	% The number of iterations to be calculated

nr_of_phaselevels=128;	% The number of phase levels that can be used
                        % OBS! If set to zero the resulting hologram will be a 
                        % pure binary amplitude hologram.

circular_beam=0;		% 1 indicates that cross section of the incoming beam is circular
                        % 0 means that the whole hologram is illuminated

circular_radie=0.9*N/2;	% The radie of the beam （if it is circular）                           


% ********** Plane 1 - The hologram plane

% Create the original electric field in plane 1 as a homogen field
% with the amplitude set to one and the phase set to zero （a plane wave
% moving along the optical axis）


% Calculate the circular cross section of the beam that is used if the variable
% “circular_beam“ is set to one

X_pos=ones（N1）*（1:N）;
Y_pos=rot90（X_pos-1）;

E_beam = ones（N）.*（（N/2-X_pos）.^2+（N/2-Y_pos）.^2
if circular_beam
    E_plane_1 = E_beam;
else
    E_plane_1 = ones（N）;
end


% ********** Plane 2 - The far field
I_plane2=zeros（N）;							% The desired intensity distribution
                                            % in plane 2

% Just a single pixel is set => The result （a plane wave） is given in only one iteration                                              
%I_hologram（N/2N/2）=1;

% It is possible to use a ordinary *.bmp picture as intensity distribution
% Just change the filename to use another image! 
% OBS! The new image has to be of the correct size!

% A binary “FOI“
%
[bild map]=imread（‘FOI.bmp‘）; 

% Check the size of the input image!
if size（bild1）~=N
    disp（[‘The input image does not have ‘ num2str（N） ‘ rows!‘]）;
    return;
end
if size（bild2）~=N
    disp（[‘The input image does not have ‘ num2str（N） ‘ columns!‘]）;
    return;
end 

I_plane_2_goal = ones（N）-double（bild）;
I_plane_2_goal_total = sum（sum（I_plane_2_goal））;

E_plane_2_goal = sqrt（I_plane_2_goal）;	% The E-field in the far field （the phase
                                        % doesn‘t matter）


% --------------------