diff --git a/Accordion_2.m b/Accordion_2.m
new file mode 100644
index 0000000..b3e21a7
--- /dev/null
+++ b/Accordion_2.m
@@ -0,0 +1,69 @@
+clear;
+
+% Elementary constants
+c = 299792458;                  %Speed of light [m/s]
+h = 6.626070040e-34;            %Planck constant [J/Hz]
+hbar = h/2/pi; 
+u0 = 4*pi*1e-7;                 %Vacuum permeability [N/A^2]
+e0 = 1/u0/(c^2);                %Vacuum permittivity [A^2*s^4/kg/m^3]
+qe = 1.6021766208e-19;          %Elementary charge [C]
+G = 6.67408e-11;                %Gravitational constant [m^3/kg/s^2]
+kB = 1.38064852e-23;            %Boltzmann constant [J/K]
+me = 9.10938356e-31;            %Electron rest mass [kg]
+mp = 1.672621898e-27;           %Proton mass [kg]
+uB = qe*hbar /2/me;             %Bohr magneton [J/T]
+uN = qe*hbar /2/mp;             %Nuclear magnton [J/T]
+alpha = u0*qe^2*c/(2*h);        %Fine structure constant
+u = 1.660539040e-27;            %Atomic mass unit [kg]
+a0 = 4*pi*e0*hbar^2/me/qe^2;    %Bohr radius [m]
+Ry = 13.605693009;              %Rydberg energy [eV]
+
+%Dysprosium properties
+mDy = 164*u;                    %Dysprosium mass [kg]
+abuDy164 = 0.2826;              %Dy164 abundance
+Polar1064 = 193*4*pi*e0*a0^3;   %Dynamical polarizability of Dy for 1064 and 532 [A^2s^4/kg]
+Polar532 = 350*4*pi*e0*a0^3;
+
+%Laser properties
+lambda = 532e-9;                %Laser wavelength [m]
+k = 2*pi/lambda;                %Wave number of laser [m^-1]
+P0 = 10;                         %Power [W]
+
+%Coordinate rotation y-axis for E1
+theta= 4/180*pi;                       %Half angle between the beams
+x = linspace(-3000,3000,1501)*1e-6;        %Definition of coordinate System
+y = 0;
+z = linspace(-300,300,1501)*1e-6;
+
+[X,Z]=meshgrid(x,z);
+xt= X.*cos(theta)+Z.*sin(theta);
+yt = y;
+zt= -X.*sin(theta)+Z.*cos(theta);
+
+%Gaussian beam propagating in xt direction
+w0 = 90e-6;                      %Beam waist at x=xt=0 (center of the trap)
+xr = pi*w0.^2/lambda;            %Rayleigh length i xt-direction
+wxt = w0.*sqrt(1+(xt./xr).^2);     %Radius to 1/e
+P = P0/2;
+
+E0 = sqrt(2*P./(pi*wxt.^2*c*e0)); %Electric field at time 0 , x=-f origin
+E1 = E0*w0./wxt.*exp(-(zt.^2+yt.^2)./wxt.^2).*exp(-1i*k.*xt); %Electric field 1
+
+%Coordinate rotation y-axis for E2
+theta= -theta;                       %Half angle between the beams
+xt= X.*cos(theta)+Z.*sin(theta);
+yt = y;
+zt= -X.*sin(theta)+Z.*cos(theta);
+
+E2 = E0*w0./wxt.*exp(-(zt.^2+yt.^2)./wxt.^2).*exp(-1i*k.*xt); %Electric field 1
+
+%Interference of the two beams
+Itot = abs(E1+E2).^2 ;             %Intensity of the interference pattern
+
+figure
+imagesc(x,z,Itot)
+set(gca,'YDir','normal')
+
+
+
+