     %%%%%%%%%%%%%   Written by Sarvenaz Memarzadeh 7/10/16  %%%%%%%%%%%%%%
clc
clear all
close all

%% Constant Parameters
q = 1.6*10^-19;
d = 0.048197;
b = 0.034274;
c0 = 0.016976;
dc = 0.001;
d1 = 0.001;
M = 9.109*10^-31;
c = 2.9979*10^8;
N = 100;
gama = 1.02;
beta = 0.195;
betaG = 0.199;
k0s = 200; 
%%
s = [-0.15:10^-5:0.15];
fs = exp(-s.^2/d^2) .* (sech(s/b) + c0 *sinh(s/b).^2);
FS = @(s) exp(-s.^2/d^2).^2 .* (sech(s/b) + c0 *sinh(s/b).^2).^2; % =fs^2
ls = integral (FS, -0.15, 0.15);
fsol = 1/k0s/ls;
fsol = fsol*100; %m.cm
a = q/(2*M*c*betaG);
B0s = sqrt(k0s) *2*M*c*gama*beta/q;
KS = @(s) exp(-s.^2/d^2).^2 .* (sech(s/b) + c0 *sinh(s/b).^2).^2*(a*B0s)^2;
ks0 = KS(0)/10^8; % 1/m^2
ks = (a*fs*B0s).^2;
HE_c = ks0 * (StepFun(s+dc/2)+ StepFun(-s+dc/2)-1);

for jj = 0:1:N;
    pha1 = StepFun (s - (dc +jj*d1)+d1/2);
    pha2 = StepFun(-1*s + dc+jj*d1+d1/2);
    pha3 = StepFun(-1*s - (dc+jj*d1)+d1/2);
    pha4 = StepFun (s + (dc +jj*d1+d1/2));
    HE_r(jj+1,:) = KS(dc + jj*d1).*(pha1 + pha2 -1); 
    HE_l(jj+1,:) = KS(dc + jj*d1).*(pha3 + pha4 -1); 
end

HE = sum(HE_c)+ sum(HE_r)+ sum(HE_l);
figure; plot(s,HE); hold on ; plot(s,ks, 'linewidth', 1.5);
xlabel('s'); legend ('HE_n(s)','\kappa(s)')


