%%%%%   EnvEqn_SmApprox1 Written by : Sarvenaz Memarzadeh 7/26/2016  %%%%%%

clc
clear all
close all

%% Constant Parameter

lsol = 0.065;
kb = 0.0015;
emit = 60*10^-6;
S = 0.32;
N = 35;
k0 = 100;
L = (N+1)*S;
del = 0.001;
M = L/del; % step size
z = linspace(0,L,M);
kh_z = kh_z1 (z);
figure; plot(z,kh_z); axis([0 L 0 110]);xlabel('z'); ylabel('\kappa_h(z)')

%% Average Focusing Function

% kave = (k0/L)*integral (@ HE, 0 , L-S/2+lsol/2)
kave = 20.31;

%% Average beam Radius from Smooth Approximation

Rave = sqrt(kb/(2*kave) + sqrt((kb/(2*kave))^2 + emit^2/kave))

%% Solution of envelop equation for UNIFORM FOCUSING and PERIODIC FOCUSING

InitCon1 = [Rave 0];
[z1,Rsm] = ode45(@myfun1 ,z ,InitCon1 );
% Beam Radius and Slope for MATCHED SOLUTION from MENV matching program :
Rm = 0.0074071; RPm = 0.00012704; 
InitCon2 = [Rm RPm];
options = odeset('RelTol', 1e-6);
[z2,Rper] = ode45(@myfun2 ,z ,InitCon2,options );
figure; plot(z1/S, Rsm(:,1),'r', 'linewidth', 1.5); hold on
plot(z2/S, Rper(:,1), 'linewidth', 1.5); axis ([0 35 0 0.01]);
xlabel('z/s'); legend('Rsm(z)', 'Rper(z)','location','southeast')

%% Solution of SINGLE PARTICLE Equation for PERIODIC FOCUSING :

InitCon3 = [0.001 -0.002];
options = odeset('RelTol',1e-6,'AbsTol',1e-8);
[z,rsp2] = ode45(@myfun3 ,z ,InitCon3,options );
figure; plot(z/S, rsp2(:,1),z2/S, Rper(:,1),'g', 'linewidth', 1.5) ; axis ([0 35 -0.01 0.01]);
xlabel('z/s'); legend('r_{sp2}(z)','R_{per}(z)','location', 'southeast')

