function NIM_evanescent (action)
% BUGS:
% Ex/Ey switch doesn't work

% dbstop if error % dbclear if error
global TestEnergy N_down N_up N_x ...
    Xup Yup Xdown Ydown ...
    N_color cmap Cycle

if nargin == 0, action = 'Init'; end;
switch action

case 'Init' % Set variables

Cycle = 1;
TestEnergy = 1; % проверка закона сохранения энергии
N_down = 100; N_up = N_down*2; N_x = N_up;
x = linspace(-2,2,N_x+1);
y_up = linspace(0,2,N_up+1);
y_down = linspace(-1,0,N_down+1);y_down = y_down(1:N_down);
[Xup,Yup] = meshgrid(x,y_up);
[Xdown,Ydown] = meshgrid(x,y_down);

q0 = 1.3;
Alpha0 = 1/2;%.33 ;%угол падения (*pi/2) в радианах
width0 = .75;%ширина пучка
n_ref0 = 0.5;%коэффициент преломления
r_down0 = 1.5; % яркость
show = [1 1 1]; % down, reflection, transmission
N_color = 512;

% Init Graph%
hL = findobj(0,'Name','Light','WindowStyle','docked');
    if isempty(hL), figure('units', 'normalized','position',[0 .035 .8 .888],...
    'Name','Light','NumberT','off',...
    'ButtonDownFcn','NIM_evanescent(''Cycle'')',...
    'CloseRequestFcn','NIM_evanescent(''Close'')');
    else figure(hL); delete(allchild(hL)); 
    end;
surf([Xdown;Xup],[Ydown;Yup],zeros(N_down+N_up+1,N_x+1),'tag','LightSurf');
shading('interp');
axis equal;
axis([min(x) max(x) min(y_down) max(y_up)]); view(2);

hL = findobj(0,'Name','Vectors','WindowStyle','docked');
    if isempty(hL),figure('units', 'normalized','position',[.8 .035 .2 .2],...
    'Name','Vectors','NumberT','off',...
    'ButtonDownFcn','NIM_evanescent(''Cycle'')',...
    'CloseRequestFcn','NIM_evanescent(''Close'')');
    else figure(hL); delete(allchild(hL));
    end;
  
    % Schema
surf([-1 1;-1 1],[-1 -1;1 1],zeros(2),cat(3,ones(2),ones(2),zeros(2)));
hold on; rotate3d on;
set(gca,'XTickLabel',[],'YTickLabel',[],'ZTickLabel',[])
alpha(.3); 
shading('interp');
axis equal;
axis([-1 1 -1 1 -1 1]); view(3);
set(gca,'Units','normalized','Position',[0 0 1 1]);

plot([0 1],[0 0],'g-','LineWidth',2,'tag','k_down');
plot([0 1],[0 0],'g-','tag','k_up');
plot([0 1],[0 0],'g-','tag','k_tr');

plot([0 1],[0 0],'m-','LineWidth',2,'tag','E_down');
plot([0 1],[0 0],'m-','tag','E_up');
plot([0 1],[0 0],'m-','tag','E_tr');

plot([0 1],[0 0],'b-','LineWidth',2,'tag','H_down');
plot([0 1],[0 0],'b-','tag','H_up');
plot([0 1],[0 0],'b-','tag','H_tr');
    
hL = findobj(0,'Name','Buttons');
    if isempty(hL),figure('units', 'normalized','position',[.8 .312 .2 .639],...
    'Name','Buttons','NumberT','off',...
    'ButtonDownFcn','NIM_evanescent(''Cycle'')',...
    'CloseRequestFcn','NIM_evanescent(''Close'')');
    else figure(hL); delete(allchild(hL));
    end;
  
% button's

%проверка закона сохраненеия энергии
    if TestEnergy,
    uicontrol('Style','text', 'String', '[r^2, t^2*n] = ',...
        'Units', 'normalized','Position',[.05 .89 .9 .09],'Tag','energ');
    end;

%угол
uicontrol('Style','slider','String','Alpha',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Units','normalized','Position',[.05 .8 .9 .03],...
    'Min',-1,'Max',1,'SliderStep',[.005 .05],...
    'Value',Alpha0,'Tag','alphaV');
uicontrol('Style','text','String',['angle = ' num2str(Alpha0) '*pi/2 рад'],...
    'Units','normalized','Position',[.05 .76 .9 .03],'Tag','alpha');

%коэффициент преломления
uicontrol('Style','slider','String','nref',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Units','normalized','Position',[.05 .71 .9 .03],...
    'Min',-5,'Max',5,'SliderStep',[.001 .1],...
    'Value',n_ref0,'Tag','nrefV');
uicontrol('Style','text','String',['nref = ' num2str(n_ref0)],...
    'Units','normalized','Position',[.05 .67 .9 .03],'Tag','nref');

%|k|
uicontrol('Style','slider','String','k_modul',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Units','normalized','Position',[.05 .62 .9 .03],...
    'Min',.1,'Max',3,'SliderStep',[.01 .1],...
    'Value',q0,'Tag','k_modulV');
uicontrol('Style','text','String',['|k| = ' num2str(q0)],...
    'Units','normalized','Position',[.05 .58 .9 .03],'Tag','k_modul');

%intensity
uicontrol('Style','slider','String','intensity',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Units','normalized','Position',[.05 .53 .9 .03],...
    'Min',1e-10,'Max',10,'SliderStep',[.01 .1],...
    'Value',r_down0,'tag','r_downV');
uicontrol('Style','text','String',['intensity = ' num2str(r_down0)],...
    'Units', 'normalized','Position',[.05 .49 .9 .03],'tag','r_down');

%ширина пучка
uicontrol('Style','slider','String','width',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Units','normalized','Position',[.05 .44 .9 .03],...
    'Min',.05,'Max',3,'SliderStep',[1/59 .2],...
    'Value',width0,'tag','widthV');
uicontrol('Style','text','String',['width = ' num2str(width0)],...
    'Units', 'normalized','Position',[.05 .4 .9 .03],'tag','width');

% поляризация
hbg = uibuttongroup(...
    'SelectionChangeFcn','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Units','normalized','Position',[.05 .26 .9 .12],'tag','polarization');
uicontrol('Style','Radio','String','No Bruster, s-pol, TM',...
    'Units','normalized','Position',[.05 5/8 .9 .25],...
    'parent',hbg,'HandleVisibility','off');
uicontrol('Style','Radio','String','Bruster, p-pol, TE',...
    'Units','normalized','Position',[.05 3/8 .9 .25],...
    'parent',hbg,'HandleVisibility','off');
uicontrol('Style','Radio','String','Natural light',...
    'Units','normalized','Position',[.05 1/8 .9 .25],...
    'parent',hbg,'HandleVisibility','off');

%waves
uicontrol('Style','checkbox','Tag','down','Units','normalized',...
    'Position',[.05 .21 .06 .025],'Value',show(1),...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;');
uicontrol('Style','text','String','Incident wave','Units','normalized',...
    'Position',[.15 .21 .8 .025]);
uicontrol('Style','checkbox','Tag','up','Units','normalized',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Position',[.05 .17 .06 .025],'Value',show(1));
uicontrol('Style', 'text', 'String', 'Reflected wave', 'Units','normalized',...
    'Position', [.15 .17 .8 .025]);
uicontrol('Style','checkbox','Tag','tr','Units','normalized',...
    'Callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'Position',[.05 .13 .06 .025],'Value',show(1));
uicontrol('Style', 'text', 'String', 'Transmitted wave', 'Units','normalized',...
    'Position', [.15 .13 .8 .025]);

%color
cmap_names = {'hot','hsv','jet','gray','bone','copper','cool'};
cmap = feval(cmap_names{1},N_color);
uicontrol('Style', 'popup', 'Tag', 'cmap',...
      'String', cmap_names, 'callback',...
      ['global cmap N_color Cycle;' ... 
      'cmap_names = get(gcbo,''String'');' ...
     'cmap = feval(cmap_names{get(gcbo,''Value'')},N_color);' ...
      'if ~Cycle, NIM_evanescent(''Cycle''); end;'], ...
      'Units', 'normalized', 'Position', [.05 .08 .9 .03]);


% прокрутка фазы, мультик
uicontrol('style','slider','Tag','phase','Value',0,...
    'Units','normalized','position',[.05 .02 .9 .02],...
    'callback','global Cycle; if ~Cycle, NIM_evanescent(''Cycle''); end;',...
    'UserData',Cycle);
NIM_evanescent ('Cycle')

    
case 'Cycle'

    Cycle = 1;
    while Cycle
        if strcmp( get(gcf,'SelectionType'),'extend')
            if get(findobj('Tag','phase'),'UserData')==3,
                delete(findobj(0,'Name','Light'));
                delete(findobj(0,'Name','Buttons'));
                delete(findobj(0,'Name','Vectors'));
                return;
            end;
            set(findobj('Tag','phase'),'UserData',...
                1-get(findobj('Tag','phase'),'UserData'));
            set(gcf,'SelectionType','normal');
        elseif strcmp( get(gcf,'SelectionType'),'open')
            set(findobj('Tag','phase'),'UserData',4); %#ok<*NODEF> % record
            set(gcf,'SelectionType','normal');
            if ~isdir('!Movie'), mkdir('!Movie'); end;
            hf = findobj(0,'Name','Light'); figure(hf);
            set(hf,...
                'PaperType','a4letter',...
                'PaperPositionMode', 'manual',...
                'PaperUnits', 'centimeters',...
                'PaperOrientation', 'portrait',...
                'PaperPosition', [1 21 6 4]);
            set(hf,'renderer','zbuffer');
            axis off;
        end;
    % 'UserData' = 1 - movie; 0 - frame by frame; 4+ - record; 3 - close;
    if get(findobj('Tag','phase'),'UserData'),
        set(findobj('Tag','phase'),'Value',...
            get(findobj('Tag','phase'),'Value')+.01);
    else
        Cycle = 0;
    end;
    UD = get(findobj('Tag','phase'),'UserData');
    if UD>=4,
        hf = findobj(0,'Name','Light');
        print (hf,'-dpng','-r150','-noui',...
            ['!Movie\MatLab'...
            num2str(UD-3,'%03d') '.png']);
        if UD >= 53,
            set(findobj('Tag','phase'),'UserData',1);
            set(hf,'renderer','OpenGL');
        else
            set(findobj('Tag','phase'),'UserData',UD+1);
        end;
    end;
    
% считывание переменных
Polarisation = ... % 1 - TM, 0 - TE (Bruster)
    strcmp(get(get(findobj('Tag','polarization'),'SelectedObject'),'String'),...
    {'No Bruster, s-pol, TM',...
    'Bruster, p-pol, TE',...
    'Natural light'});
q = get(findobj('Tag','k_modulV'), 'Value');%|k|падающей волны
set(findobj('Tag','k_modul'),'String',...
    ['|k| = ' num2str(q,2)]);
Alpha = pi/2*get(findobj('Tag','alphaV'),'Value');%угол падения в радианах
set(findobj('Tag','alpha'),'String',...
    ['angle = ' num2str(Alpha*2/pi,2) '*pi/2 рад']);
width = get(findobj('Tag','widthV'), 'Value');%ширина пучка
set(findobj('Tag','width'),'String',...
    ['width =' num2str(width,2)]);
n_ref = get(findobj('Tag','nrefV'),'Value');% коэффициент преломления
set(findobj('Tag','nref'),'String',...
    ['nref = ' num2str(n_ref,2)]);
    % if n_ref == 0, n_ref = 1e-5; end;
r_down = get(findobj('Tag','r_downV'),'Value'); % яркость
set(findobj('Tag','r_down'),'String',...
    ['intensity = ' num2str(r_down,2)]);
show(1) = get(findobj('Tag','down'),'Value');
show(2) = get(findobj('Tag','up'),'Value');
show(3) = get(findobj('Tag','tr'),'Value');
hs = findobj('Tag','phase');
phi = get(hs,'Value');
    if phi >= 1, set(hs,'Value',eps); end;
    if phi <= 0, set(hs,'Value',1); end;

% расчет переменных
Psi = n_ref/abs(n_ref)*asin(sin(Alpha)/(n_ref)); % угол преломления
q_ref = q*n_ref;%|k| преломленной волны

k = q*[sin(Alpha)  -cos(Alpha)];%к вектор падающей волны
k_ref = [k(1)  -n_ref/abs(n_ref)*sqrt(q_ref^2-k(1)^2)]; %к вектор преломленной волны
n = abs(n_ref);

r_up = 0;
r_tr = 0;
if abs(sin(Alpha))<n,
    if Polarisation(1), % No Bruster, s-pol, TM
        r_up = r_down*(cos(Alpha)-n*cos(Psi))/(cos(Alpha)+n*cos(Psi));
        %коэффициент Френеля для отраженной волны
        r_tr = r_down*2*cos(Alpha)/(cos(Alpha)+n*cos(Psi));
        %коэффициент Френеля для прошедшей
    elseif Polarisation(2) % Bruster, p-pol, TE
        r_up = r_down*(n*cos(Alpha)-cos(Psi))/(n*cos(Alpha)+cos(Psi));
        r_tr = r_down*2*cos(Alpha)/(n*cos(Alpha)+cos(Psi));
    elseif Polarisation(3) % Natural light
        r_up = r_down/sqrt(2)*sqrt(...
            (cos(Alpha)-n*cos(Psi))^2/(cos(Alpha)+n*cos(Psi))^2+...
            (n*cos(Alpha)-cos(Psi))^2/(n*cos(Alpha)+cos(Psi))^2);
        r_tr = r_down*sqrt(2)*sqrt(cos(Alpha)^2/(cos(Alpha)+n*cos(Psi))^2+...
            cos(Alpha)^2/(n*cos(Alpha)+cos(Psi))^2);
    end;
    width_ref = width * cos(Psi)/cos(Alpha);% ширина преломленного луча
    E_tr = r_tr*exp(-((Xdown*k_ref(2)-Ydown*k_ref(1))/q_ref/width_ref).^2)...
        .*real(exp(1i*2*pi*(Xdown*k_ref(1)+Ydown*k_ref(2)-phi))); %прошедшая волна
    D = 0; % Goos-Hanchen
else % полное внутреннее отражение
    if Polarisation(1)||Polarisation(1), % No Bruster, s-pol, TM !!! cheat natural light
        D = sign(n_ref)*1/pi/q*... % Goos-Hanchen
            sin(Alpha)*cos(Alpha)^2/(1-n_ref^2)...
            /sqrt(sin(Alpha)^2-n_ref^2);
    elseif Polarisation(2) % Bruster, p-pol, TE
        D = sign(n_ref)*1/pi/q*...
            n_ref^2*sin(Alpha)*cos(Alpha)^2 ...
            /(n_ref^2*cos(Alpha)^2+sin(Alpha)^2-n_ref^2) ...
            /sqrt(sin(Alpha)^2-n_ref^2);
    end;
    Psi = 0;
    r_up = r_down; r_tr = 2*r_down;%/n_ref; 
    width_ref = width/cos(Alpha);% ширина преломленного луча
    % Psi = NaN;%n_ref/abs(n_ref)*asin(sin(Alpha)/(n_ref)); % угол преломления
    k_ref(2) =  -sqrt(q_ref^2-k_ref(1)^2); %к вектор преломленной волны поперек границы
    E_tr = r_tr*real(exp(-(max(0,abs(Xdown)-abs(D)/2/cos(Alpha))/width_ref).^2)...
        .*real(exp(1i*2*pi*(Xdown*k_ref(1)+Ydown*k_ref(2)-phi)))); %прошедшая волна
end;
% phase_trans = 1/2;
E_down = r_down*exp(-((Xup*k(2)-Yup*k(1)-D/2)/q/width).^2)...
    .*cos(2*pi*(Xup*k(1)+Yup*k(2)-phi)); %падающаю волна
E_up = r_up*exp(-((Xup*k(2)+Yup*k(1)+D/2)/q/width).^2)...
    .*cos(2*pi*(Xup*k(1)-Yup*k(2)-phi));%отраженная волна

if ~show(1), E_down = 0*E_down; end;
if ~show(2), E_up = 0*E_up; end;
if ~show(3), E_tr = 0*E_tr; end;
C_int = [n*E_tr.^2;(E_up+E_down).^2];

CInd = floor(1+atan(C_int*pi/2)*2/pi*(N_color-1));
CDat = reshape(cmap(CInd,:),...
    [N_up+N_down+1,N_x+1,3]);

set(findobj('Tag','LightSurf'),'ZData',C_int);
set(findobj('Tag','LightSurf'),'CData',CDat,'CDataMapping','direct');

set(findobj('tag','k_down'),'XData',[0 -sin(Alpha)],'ZData',[0 cos(Alpha)]);
set(findobj('tag','k_up'),'XData',[0 sin(Alpha)],'ZData',[0 cos(Alpha)]);
set(findobj('tag','k_tr'),'XData',[0 sin(Psi)],'ZData',[0 -cos(Psi)]);

if Polarisation(1), % No Bruster, s-pol, TM
    set(findobj('tag','E_down'),'XData',[0 0],'ZData',[0 0],'YData',[0 -1]);
    set(findobj('tag','E_up'),'XData',[0 0],'ZData',[0 0],'YData',[0 -r_up/r_down]);
    set(findobj('tag','E_tr'),'XData',[0 0],'ZData',[0 0],'YData',[0 -r_tr/r_down]);

    set(findobj('tag','H_down'),'XData',[0 -cos(Alpha)],...
        'YData',[0 0],'ZData',[0 -sin(Alpha)]);
    set(findobj('tag','H_up'),'XData',[0 r_up/r_down*cos(Alpha)],...
        'YData',[0 0],'ZData',[0 -r_up/r_down*sin(Alpha)]);
    set(findobj('tag','H_tr'),'XData',[0 -n_ref*r_tr/r_down*cos(Psi)],...
        'YData',[0 0],'ZData',[0 -n_ref*r_tr/r_down*sin(Psi)]);
elseif Polarisation(2) % Bruster, p-pol, TE
    set(findobj('tag','H_down'),'XData',[0 0],'ZData',[0 0],'YData',[0 -1]);
    set(findobj('tag','H_up'),'XData',[0 0],'ZData',[0 0],'YData',[0 -r_up/r_down]);
    set(findobj('tag','H_tr'),'XData',[0 0],'ZData',[0 0],'YData',[0 -r_tr/r_down*n_ref]);

    set(findobj('tag','E_down'),'XData',[0 cos(Alpha)],...
        'YData',[0 0],'ZData',[0 sin(Alpha)]);
    set(findobj('tag','E_up'),'XData',[0 -r_up/r_down*cos(Alpha)],...
        'YData',[0 0],'ZData',[0 r_up/r_down*sin(Alpha)]);
    set(findobj('tag','E_tr'),'XData',[0 r_tr/r_down*cos(Psi)],...
        'YData',[0 0],'ZData',[0 r_tr/r_down*sin(Psi)]);
elseif Polarisation(3) % Natural
    set(findobj('tag','E_down'),'XData',[0 0],'ZData',[0 0],'YData',[0 0]);
    set(findobj('tag','E_up'),'XData',[0 0],'ZData',[0 0],'YData',[0 0]);
    set(findobj('tag','E_tr'),'XData',[0 0],'ZData',[0 0],'YData',[0 0]);
    set(findobj('tag','H_down'),'XData',[0 0],'ZData',[0 0],'YData',[0 0]);
    set(findobj('tag','H_up'),'XData',[0 0],'ZData',[0 0],'YData',[0 0]);
    set(findobj('tag','H_tr'),'XData',[0 0],'ZData',[0 0],'YData',[0 0]);
end;

if TestEnergy,
    Transition = n*cos(Psi)/cos(Alpha)*r_tr^2/r_down^2;
    if isnan(Transition), Transition = 0; end;
    set(findobj('Tag','energ'),'string',...
        ['[r, t] = [' num2str(r_up/r_down,3) ', ' num2str(r_tr/r_down,3) ']' ...
        char([10 13]) '[r^2, t^2*n] = ['...
        num2str(r_up^2/r_down^2,3) ', ' num2str(Transition,3) ']' ...
        char([10 13]) 'E/E0 = ' num2str(r_up^2/r_down^2+Transition,3)]);
else
    if abs(1-(r_up^2+n*cos(Psi)/cos(Alpha)*r_tr^2)/r_down^2)>1e-10,
        warning('Enegry is not conserved!'); %#ok<WNTAG>
    end;
end;

% pause(.02); 
drawnow;
    end;
case 'Close'
    if get(findobj('Tag','phase'),'UserData')==1,
        set(gcf,'SelectionType','extend');
        set(findobj('Tag','phase'),'UserData',3);
    else
        delete(findobj(0,'Name','Light'));
        delete(findobj(0,'Name','Buttons'));
        delete(findobj(0,'Name','Vectors'));
    end;
end;
end
