help with curve fitting sinc^2(x)
    12 views (last 30 days)
  
       Show older comments
    
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
      
      
 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.
Answers (0)
See Also
Categories
				Find more on Interpolation in Help Center and File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
