Notes

CORONASIMULATE Simulate coronagraph and Gerchberg-Saxton algorithm
%
% A simulation of a coronagraph and the Gerchberg-Saxton algorithm, in the
% context of NASA's Roman Space Telescope, developed to help teach ENCMP
% 100 Computer Programming for Engineers at the University of Alberta. The
% program saves data to a MAT file for subsequent evaluation.

%{
Copyright (c) 2021, University of Alberta
Electrical and Computer Engineering
All rights reserved.

Student name:Abdullah
Student CCID:Sherwani
Others:

To avoid plagiarism, list the names of persons, other than a lecture
instructor, whose code, words, ideas, or data were used. To avoid
cheating, list the names of persons, other than a lab instructor, who
provided substantial editorial or compositional assistance.

After each name, including the student's, enter in parentheses an
estimate of the person's contributions in percent. Without these
numbers, adding to 100%, follow-up questions may be asked.

For anonymous sources, enter code names in uppercase, e.g., SAURON,
followed by percentages as above. Email a list of online sources
(links suffice) to a teaching assistant before submission.
%}
clear
rng('default')

im = loadImage;
[im,Dphi,mask] = opticalSystem(im,300);
[images,error] = gerchbergSaxton(im,20,Dphi,mask);

frames = getFrames(images,error);
save frames frames



% im = loadImage returns an indexed image, im, of an exoplanet (2M1207b)
% orbiting a distant star. The image, a 2D matrix, is downloaded from a
% NASA website. It should be displayed with a grayscale colour map.
function im = loadImage
path = 'https://exoplanets.nasa.gov/system/resources/';
file = 'detail_files/300_26a_big-vlt-s.jpg';
im = imread(strcat(path,file));
im = rgb2gray(im);
end

function [im,Dphi,mask] = opticalSystem(im,width)
[im,mask] = occultCircle(im,width);
[IMa,IMp] = dft2(im);
imR = uint8(randi([0 255],size(im)));
[~,Dphi] = dft2(imR);
im = idft2(IMa,IMp-Dphi);
end

function [im,mask] = occultCircle(im,width)
%this locates the center of the circle
o = size(im)/2;
%places in circle
im=insertShape(im,'filledcircle',[o(2) o(1) width/2],'color','black','opacity',1);
%makes image gray
im=rgb2gray(im);
for s = 1:2*o(2)
for t = 1:2*o(1)
if 0 == im(t,s)
%white pixel for true
mask(t,s) = 255;
else
%black pixel for false
mask(t,s) = 0;
end
end
end
end

% [IMa,IMp] = dft2(im) returns the amplitude, IMa, and phase, IMp, of the
% 2D discrete Fourier transform of an indexed image, im. The image, a 2D
% matrix, is expected to be of class 'uint8'. The phase is in degrees.
function [IMa,IMp] = dft2(im)
IM = fft2(im);
IMa = abs(IM);
IMp = rad2deg(angle(IM));
end

% im = idft2(IMa,IMp) returns an indexed image, im, of class 'uint8' that
% is the inverse 2D discrete Fourier transform (DFT) of a 2D DFT specified
% by its amplitude, IMa, and phase, IMp. The phase must be in degrees.
function im = idft2(IMa,IMp)
IM = IMa.*complex(cosd(IMp),sind(IMp));
im = uint8(ifft2(IM,'symmetric'));
end

function [image,errors] = gerchbergSaxton(im,maxIters,Dphi,mask)
[IMa,IMp] = dft2(im);
image = cell(maxIters+1,1);
errors = zeros(maxIters+1,1);
for p = 0:maxIters
fprintf('Iteration %d of %dn',p,maxIters)
im = idft2(IMa,IMp+((p/maxIters)*Dphi));
image{p+1} = im;
o=size(mask);
point = 1;
%find mask pixels
for n = 1:o(2)
for i = 1:o(1)
if mask(i,n) == 255
imm(point) = im(i,n);
point = point + 1;
else
imm(point) = 0;
point = point + 1;
end
end
end
%stores then saves the values
immp = imm;
immp = im2double(immp);
immp1 = immp';
%saves error values
errors(p+1) = immp*immp1;
end
end

function frames = getFrames(images,error)
NumberFrames = numel(images);
for S = 1:NumberFrames
hold on
image([0 NumberFrames],[max(error) 0], images{S});
plot(0:S-1,error(1:S),'r-');
title("Coronagraph Simulation");
xlabel("X pixel location")
ylabel("Y pixel location")
%This Labels the title and the axes
colormap(gray);%this makes the image gray
frames(S) = getframe(gcf);
end

close all
end