help with curve fitting sinc^2(x)

22 visualizaciones (últimos 30 días)
Patrick Bevington
Patrick Bevington el 19 de Mzo. de 2016
Comentada: Star Strider el 19 de Mzo. de 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 comentario
Star Strider
Star Strider el 19 de Mzo. de 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.

Iniciar sesión para comentar.

Respuestas (0)

Categorías

Más información sobre Linear and Nonlinear Regression en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by