help with curve fitting sinc^2(x)

32 visualizzazioni (ultimi 30 giorni)
Patrick Bevington
Patrick Bevington il 19 Mar 2016
Commentato: Star Strider il 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 Commento
Star Strider
Star Strider il 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.

Accedi per commentare.

Risposte (0)

Categorie

Scopri di più su Get Started with Curve Fitting Toolbox in Help Center e File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by