How to properly scale the lorentzian curve to the data

clear all
clc
close all
RMS=[];
R=[];
P_threshold=[];
L=5; %[m] Fiber length
V1=[1];
P0=4.7960;
F1=0;
count=0;
count1=0;

%% general constants
n=1.45; %Index of refraction
eps0=8.854e-12; % [F/m] Vacuum permittivity
mu0 = 4*pi*1e-7;%[H/m] Vacuum permeability
c=2.9979e8; % [m/sec] Speed of light
Z0=sqrt(mu0/eps0);   %[Ohm] Vacuum impedance
dt = 6e-12; dz=dt*c/n;    %Spacial and Temporal step sizes.
N=round(L/dz); % Fiber length discretization
z=0:dz:(dz*N-dz);
T=20*2*L*n/c;   %time taken for 10 round trips
Nt=round(T/dt);

%% material characteristics
lambdaL=1.064e-6;            % [m] pump wavelength
WL=2*pi*c/lambdaL;       %[rad/sec]  pump frequency
GAMMAb=2*pi/17.5e-9; %[1/t] phonons decay rate(zeringue)
A=7.85e-11;   %[m^2] fiber's effective area
Va=5960;    %[m/s] Velocity of sound
Omega=2*n*Va*WL/c; %[rad/sec] (zeringue)  
WS = WL-Omega;
gammaE=1.95;  %electrostrictive coefficient
kT=300*1.3806504e-23;       %[Joule]        at room temprature. 
roh0=2201;       % [kg/m^3]  mean density for SiO2, bulk (http://www.virginiasemi.com/pdf/generalpropertiesSi62002.pdf) . 
Q=2*kT*roh0*GAMMAb/(Va^2*A);     %Noise std.
g0=(gammaE^2*WL^2)/n/c^3/Va/roh0/GAMMAb; %[m/W] (Jenkins)
sigma=(WL*gammaE)/2/n/roh0/c;    %(derivation)
q=2*WL*n/c;                       %derivation
I1_0=P0/A;  %[W/m^2] Pump intensity
LAMBDA=gammaE*Omega/c/Z0/2/Va^2;    %derivation
kappa=g0*GAMMAb*n/2/Z0/LAMBDA;      %derivation
Nz=N;

%% Coefficients
G=g0*I1_0*L;                         %Gain factor
Gth=21; %Agrawal, pp.361
Elt=repelem(0,Nz);
Elt(1)=sqrt(P0*Z0/A/2/n);
Est=repelem(0,Nz);

%%
t = (-Nt/2:1:Nt/2-1)*dt;
nu = (-Nt/2:1:Nt/2-1)*1/T;
vpi = 3;
RL=50;
V=1;
fsine=2e9;
sine = V*sin(2*pi*fsine*t);

%% White noise modulation
np = 22; % power of RF White noise signal in dBm
np_l = 10^(np/10)*1e-3; % power of RF White noise signal in Watts
nbw = 0.5e9; % bandiwdth in Hz of the RF WNS

rng('default');
rng(2);
n1 = randn(1, Nt); % generate a white gaussian noise in time domain which is delta correlated in time domain
s=nbw*sinc(pi*t*nbw);
n_LPF_T=conv(n1,s,'same');
n_lpf = sqrt(np_l*T*RL)*(n_LPF_T).*(1./sqrt(trapz(t, abs(n_LPF_T).^2))); %normalize the noise volatge signal

%%
phi=n_lpf;
ES_0t=sqrt(I1_0/2/n/c/eps0)*exp((1i*phi));  % Original signal in time 
Power=trapz(t,2*n*c*eps0*A*abs(ES_0t).^2)/T;        % Area under the curve in time domain 
FFt_EL0t=fftshift(abs(fft(ES_0t)));             % fourier transform of the original signal 
Power_FFt=T*trapz(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2);   % Area under the curve in frequency domain 
% figure;plot(nu,2*n*c*eps0*A*(FFt_EL0t/Nt).^2)

%% ---------------------------------
%% Fitting with Lorentzian model
%% ---------------------------------

y   = 2*n*c*eps0*A*(FFt_EL0t/Nt).^2;
[up, ~] = envelope(y, round(0.01*size(nu/max(nu), 2)), 'peak');

format long g
x   = nu;
xm  = max(nu)

xm = 
          83332316172.4138

x_m=std(x);
x_axis=(x)/xm;
% b1= sqrt(1/(max(y)/min(y)-1));
% b1=(3343175379.31034/xm)^2;
% b1=1/((max(x)-min(x))/10/xm)^2;
% b1=(abs(x(dsearchn(y',max(y)/10)))/xm)^2;
c1=x(dsearchn(y',max(y)));
% b1=((xm-c)^2/((max(y)/min(y))-1));
a1=max(y);
b1=abs(x(dsearchn(y',max(y)/2)));
tol=0.0001;
y_lor=a1./(((x-c1)/b1).^2+1);
% Fitting model
 fo  = fitoptions('Method', 'NonlinearLeastSquares',...
              'Lower', [tol*a1, tol*b1,tol-c1 ],...
                'Upper', [ a1, b1,  tol+c1]);
% fo  = fitoptions('Method', 'NonlinearLeastSquares' );
ft  = fittype('(a)./(((x-c)/b).^2+1)', ... 
    'independent', {'x'}, ...
    'coefficients', {'a','b', 'c'}, ...
    'options', fo);

%% Fit curve to data
[yfit, gof] = fit(x', y', ft, 'StartPoint', [a1, b1,c1])

yfit = 
     General model:
     yfit(x) = (a)./(((x-c)/b).^2+1)
     Coefficients (with 95% confidence bounds):
       a =    0.001715  (0.001703, 0.001727)
       b =   1.148e+09  (1.136e+09, 1.159e+09)
       c =      0.0001  (fixed at bound)

gof = struct with fields:
           sse: 0.0053242226215525
       rsquare: 0.483756138064637
           dfe: 161222
    adjrsquare: 0.483752935996297
          rmse: 0.000181725532855836

figure
plot(yfit, x, y), grid on, grid minor, xlim([-12e9, 12e9])
xlabel('\nu'), ylabel('y')

How do I increase the amplitude of the curve to properly match the data...