Clear Filters
Clear Filters

help with curve fitting sinc^2(x)

22 views (last 30 days)
Patrick Bevington
Patrick Bevington on 19 Mar 2016
Commented: Star Strider on 19 Mar 2016
Hi, I have written some code to plot the intensity of light seen by a diffraction grating and I want to fit the curve to some data points I have but I can't get it work. Can anyone see where might be going wrong? What ever values I pick for starting point I can't get the data to plot. The code named 'Intensity.m' is for plotting the curve I want i.e. the distribution I want, and the code named 'CurveFit.m' is to try and fit the curve to my data points.
'Intensity.m'
clear all, clc
a = 0.9*20.4e-6; % Width - amplitude, increase a -> amp smaller
d = 2.1*20.4e-6; % Seperation - peak on x, increase d -> peaks closer
l = 780e-9; % Wavelength
s = 1; % Distance
thetamax = pi/50; % Angle
theta = -thetamax : 1e-5 :thetamax;
x=s*tan(theta); % Order seperation
alpha=pi*d*sin(theta)/l;
y1=cos(alpha).^2; % Interference
beta = pi*a*sin(theta)/l;
y2 = (sin(beta)./beta).^2; % Diffraction envolpe
y = y1 .* y2; % Interference Pattern
figure(1)
plot(x,y)
title('Grating pattern');
xlabel('Distance');
ylabel('Intensity');
% axis([-0.02,0.02,0,inf])
Here I used the symbolic function in matlab to generate an equation for y, give to be
y = (l^2*cos((pi*d*sin(theta))/l)^2*sin((pi*a*sin(theta))/l)^2)/(a^2*pi^2*sin(theta)^2)
This is my attempt at the curve fitting. I need it so my data points y1,y2,y3... fit on top of the central peak and the two either side.
'CurveFit.m'
clear all, clc
y1 = [0.2685, 1, 0.2981]; % Data to fit
y2 = [0.3339, 1, 0.3290];
y3 = [0.3039, 1, 0.3012];
y4 = [0.4269, 1, 0.4354];
y5 = [0.3232, 1, 0.3608];
y6 = [0.3582, 1, 0.4291];
y7 = [0.3767, 1, 0.4491];
y8 = [0.3186, 1, 0.4245];
x = linspace(-pi,pi,3)';
aa = 20.4e-6; % Width a
bb = 20.4e-6; % Seperation d
startPoints = [aa bb];
% y = (l^2*cos((pi*d*sin(theta))/l)^2*sin((pi*a*sin(theta))/l)^2)/(a^2*pi^2*sin(theta)^2)
FitEquation = fittype('(780e-9^2 * cos( (pi * b * sin(x) )/780e-9 )^2 * sin( (pi * a * sin(x)) /780e-9 )^2 )/(a^2 * pi^2 * sin(x)^2)');
f1 = fit(x,y1,FitEquation,'Start',startPoints);
figure(99)
plot(f1,x,y1)
  1 Comment
Star Strider
Star Strider on 19 Mar 2016
Your objective function may be wrong. I did the fit with this:
% MAPPING: B(1) = a, B(2) = b
diffrax = @(B,x) (780e-9^2 * cos( (pi * B(2) .* sin(x) )/780e-9 ).^2 .* sin( (pi * B(1) .* sin(x)) /780e-9 ).^2 )./(B(1).^2 * pi^2 * sin(x).^2);
y1 = [0.2685, 1, 0.2981]; % Data to fit
y2 = [0.3339, 1, 0.3290];
y3 = [0.3039, 1, 0.3012];
y4 = [0.4269, 1, 0.4354];
y5 = [0.3232, 1, 0.3608];
y6 = [0.3582, 1, 0.4291];
y7 = [0.3767, 1, 0.4491];
y8 = [0.3186, 1, 0.4245];
D = [y1; y2; y3; y4; y5; y6; y7; y8];
D = sortrows(D, 1);
SSECF = @(B) sum((D(:,3) - diffrax(B,D(:,1))).^2);
B0 = [2E-5; 2E-5]*0.001;
[EstB, SSE] = fminsearch(SSECF, B0);
figure(1)
plot(D(:,1), D(:,3), 'bp')
hold on
plot(D(:,1), diffrax(EstB,D(:,1)), '-r')
hold off
grid
and did not get an acceptable result.
I’m not submitting this as an Answer because it isn’t one.

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!