diff --git a/Matrix.m b/Matrix.m
new file mode 100644
index 0000000..e34aef8
--- /dev/null
+++ b/Matrix.m
@@ -0,0 +1,9 @@
+function m = Matrix(Phi,Z)
+Limit=length(Phi);
+M1=[1 0;0 1];
+for j=2:Limit
+    M{j} = [cosd(Phi(j)) (1i*sind(Phi(j)))/Z(j) ; 1i*Z(j)*sind(Phi(j)) cosd(Phi(j))];
+    M1=M1*M{j};
+end
+m=M1;
+end
diff --git a/MultiLayerFilm.m b/MultiLayerFilm.m
new file mode 100644
index 0000000..3f158b8
--- /dev/null
+++ b/MultiLayerFilm.m
@@ -0,0 +1,18 @@
+function [Incident,RS,RP,TS,TP] = MultiLayerFilm(n,d,Incident,Lambda)
+BOUNDRY = length(n)-1;
+INC=Incident;
+for INCIDENCE=Incident
+%% Calculation of angles on boundries starting from first
+Theta=SnellsLaw(n,INCIDENCE);
+%% 
+Phi(2:BOUNDRY) = n(2:BOUNDRY).*d(1:BOUNDRY-1).*(2.*pi./(Lambda));
+Z_s = (2.6544e-3).*n(1:BOUNDRY+1).*cosd(Theta(1:BOUNDRY+1));
+Z_p = (2.6544e-3).*n(1:BOUNDRY+1)./cosd(Theta(1:BOUNDRY+1));
+m1=Matrix(Phi,Z_s);
+[R_s(INCIDENCE+1),T_s(INCIDENCE+1)]=R_T(m1,Z_s(1),Z_s(BOUNDRY+1));
+m2=Matrix(Phi,Z_p);
+[R_p(INCIDENCE+1),T_p(INCIDENCE+1)]=R_T(m2,Z_p(1),Z_p(BOUNDRY+1));
+end
+RS=R_s;,RP=R_p;,TS=T_s;,TP=T_p;
+rtplot(Incident,R_s,R_p,T_s,T_p);
+end
diff --git a/R_T.m b/R_T.m
new file mode 100644
index 0000000..5f0e62d
--- /dev/null
+++ b/R_T.m
@@ -0,0 +1,8 @@
+function [R,T] = R_T(m,Z_o,Z_s)
+r=(Z_o*m(1,1)+Z_o*Z_s*m(1,2)-m(2,1)-Z_s*m(2,2))/...
+            (Z_o*m(1,1)+Z_o*Z_s*m(1,2)+m(2,1)+Z_s*m(2,2));
+t=(2*Z_o)/...
+            (Z_o*m(1,1)+Z_o*Z_s*m(1,2)+m(2,1)+Z_s*m(2,2));
+T=abs(t.*t');
+R=abs(r.*r');
+end
\ No newline at end of file
diff --git a/SnellsLaw.m b/SnellsLaw.m
new file mode 100644
index 0000000..970f126
--- /dev/null
+++ b/SnellsLaw.m
@@ -0,0 +1,8 @@
+function angles = SnellsLaw(n, incident)
+Theta=zeros(1,length(n));
+Theta(1)=incident;
+for j=2:length(n)
+    Theta(j) = asind((n(j-1)/n(j))*sind(Theta(j-1)));
+end
+angles=Theta;
+end
\ No newline at end of file
diff --git a/norm2unity.m b/norm2unity.m
new file mode 100644
index 0000000..7a84957
--- /dev/null
+++ b/norm2unity.m
@@ -0,0 +1,3 @@
+function norm = norm2unity(input)
+norm=(input-min(input))/(max(input)-min(input))
+end
diff --git a/rtplot.m b/rtplot.m
new file mode 100644
index 0000000..7c10558
--- /dev/null
+++ b/rtplot.m
@@ -0,0 +1,13 @@
+function Output = rtplot(Incident,R_s,R_p,T_s,T_p)
+Incident=Incident+1;
+plot(Incident,T_s,'r');
+hold on
+plot(Incident,R_s,'r');
+hold on
+plot(Incident,T_p,'b');
+plot(Incident,abs(R_p),'b');
+hold off
+%[K,I] = min(R_p);
+%X_b=I.*(ones(1,length(Incident)));
+%plot(X_b,R_p,'--');
+end
\ No newline at end of file