function [z,t,theta,phi,u,intensity] = lc3kdt(volts)
%LC3K Calculates liquid crystal orientation in a 1D cell.
%   [Z,T,THETA,PHI,U,INTENSITY] = LC3KDT(VOLTS) returns the 
%   tilt (THETA), twist (PHI) and potential (U) as a 
%   function of Z [um] and time T [ms] for 1D cell with a 
%   potential difference of VOLTS applied across it.
%
%   This is a 3 elastic constant formulation, using a polar
%   representation of the LC orientation. For large voltages
%   the twist becomes degenerate (Jacobian matrix becomes 
%   singular).
%
%   Supports alignment layers of finite thickness. The left
%   and right alignment layers may be of different 
%   thicknesses and relative permittivities. The thickness
%   may also be set to zero.
%
%   Anchoring may be strong or weak at the left and right
%   alignment layers. Set W = Inf for strong anchoring.
%   The anchoring energy is assumed isotropic and takes 
%   the form Fs = -W/2*(n.n0)^2.
%
%   A simple Jones code calulates the optical transmittance, 
%   returned in INTENSITY. Optical properties are defined in
%   the function TRANSMITTANCE.
%
%   Copyright (C) 2010 Richard James (rwjames@gmail.com)
%
%   See also PDEPE, PDEX4.

    max_time = 10e-3; % maximum simulation time [s]
   
    K11 = 12.5e-12; % elastic constant [N] (splay)
    K22 =  6.5e-12; % elastic constant [N] (twist)
    K33 = 16.5e-12; % elastic constant [N] (bend)
    q0 = 0;         % chirality 2*pi/pitch
    epslow = 5.2;   % epsilon perpendicular LC
    deleps = 13.8;  % dielectric anisotropy LC
    pretiltl = 6*pi/180;     % pretilt at left alignment layer
    pretiltr = -5*pi/180;     % pretilt at right alignment layer
    totaltwist = 90*pi/180; % diff. in twist between the two alignment layers
    e0 = 8.85418782e-12; %permitivity of free space

    d_alignl = 0.05;   % left alignment layer thickness  [microns] can be 0
    d_alignr = 0.05;   % right alignment layer thickness [microns] can be 0
    d_LC = 1.0;%-d_alignl-d_alignr;   % liquid crystal layer thickness [microns]

    eps_alignl = 3;    % left aligment layer relative permittivity
    eps_alignr = 4;    % right aligment layer relative permittivity
    
    % left and right anchoring strengths
    Wl = 500e-6; % Weak = ? [J/m^2] or Strong = Inf
    Wr = 2000e-6; % Weak = ? [J/m^2] or Strong = Inf
    
    gamma1 = 0.1232;

    z = linspace(0, d_LC, 41);
%	t = linspace(0, max_time/d^2, 51);
    t = linspace(0, max_time/1e-6^2, 51);
    
    u0 = getu0(z); % get initial electric potential distribution
    m = 0;     
    sol = pdepe(m,@pde_3k,@ic_fun,@bc_fun,z,t);
   
%   z = z.*d*1e6;   % [um] undimensionalise
%   t = t.*d^2*1e3; % [ms]
    t = t.*1e-6^2*1e3;
    [zz,tt] = meshgrid(z, t);
    theta = sol(:,:,1);
    phi = sol(:,:,2);
    u = sol(:,:,3);

    subplot(2,2,1),
    surf(zz,tt,theta*180/pi);
    xlabel('z [\mum]'); ylabel('time [ms]')
    zlabel('tilt [degrees]')

    subplot(2,2,2),
    surf(zz,tt,phi*180/pi);
    xlabel('z [\mum]'); ylabel('time [ms]')
    zlabel('twist [degrees]')

    subplot(2,2,3),
    surf(zz,tt,u);
    xlabel('z [\mum]'); ylabel('time [ms]')
    zlabel('potential [V]')
   
    % calculate transmittance vs. time
    intensity = zeros(1,size(theta,1));
    for it=1:size(theta,1)
    	intensity(it) =  transmittance(z, theta(it,:), phi(it,:));
    end
    subplot(2,2,4),
    plot(t, intensity, 'ro-')
    xlabel('time [ms]')
    ylabel('transmittance')

    % Oseen Frank free energy and Laplace equation
    % ------------------------------------------------------------
    function [c,f,s] = pde_3k(z,t,u,DuDx)%,params)
        % c du/dt = d/dz [ f(u,du/dz) ] + s(u,du/dz)
        costheta1 = cos(u(1));
        sintheta1 = sin(u(1));
        theta2 = DuDx(1); % the spatial derivative theta'
        phi2 = DuDx(2); % the spatial derivative phi'
        u2 = DuDx(3); % the spatial derivative u'

        c = [gamma1; gamma1*costheta1^2; 0];

        % maple generated (see lc3k_dt.mws) 
        % split Euler-Lagrange into two contributions:
        % dF/d(theta) (only tilt contibutes)
        s(1) = -K33*theta2^2*costheta1*sintheta1+K11*costheta1*theta2^2* ...
            sintheta1-2*K33*costheta1^3*phi2^2*sintheta1+K33*costheta1* ...
            phi2^2*sintheta1-2*K22*costheta1*phi2*q0*sintheta1+2*K22* ...
            costheta1^3*phi2^2*sintheta1+e0*deleps*u2^2*costheta1*sintheta1;
        s(2) = 0;
        s(3) = 0;

        % derivative terms: dF/d(d theta/dz) & dF/d(d phi/dz)
        f(1) = -K33*theta2*costheta1^2+K11*costheta1^2*theta2+K33*theta2;
        f(2) = -K33*costheta1^4*phi2-K22*costheta1^2*q0 ...
            +K22*costheta1^4*phi2+K33*costheta1^2*phi2;
        f(3) = e0*(epslow + deleps*sintheta1^2)*u2;
    end % function pde_3k

    % boundary conditions 
    % ------------------------------------------------------------
    function [pl,ql,pr,qr] = bc_fun(zl,ul,zr,ur,t)
        [pl,ql] = surface_fun(-Wl, ul(1), ul(2), pretiltl, 0, ul(3), 0);
        [pr,qr] = surface_fun( Wr, ur(1), ur(2), pretiltr, totaltwist, ur(3), volts);

        % Neumann boundary conditions for alignment layers
        % note that for the Neumann boundary condition cases, the natural
        % boundary condition implied by the functional is satified, if we set
        % the appropriate ql and qr value to \pm1. See PDEX4 for example.
        if d_alignl~=0
            pl(3) = e0*eps_alignl*(0-ul(3))/d_alignl;
            ql(3) = 1;
        end
        if d_alignr~=0
            pr(3) = e0*eps_alignr*(volts-ur(3))/d_alignr;
            qr(3) = -1;
        end
    end % function bc_fun

    function [p, q] = surface_fun(W, theta1, phi1, theta0, phi0, u1, u0)
        if isinf(W), % strong anchoring
            p = [theta1-theta0; phi1-phi0; u1-u0];
            q = [0; 0; 0]; % derivative
        else % weak anchoring
            % anchoring energy is of the form -W/2*(n.n0)^2, change these 
            % terms if you want anisotropic anchoring
            tmp = -W*1e-6*(cos(theta1)*cos(theta0)*(cos(phi1)*cos(phi0)+ ...
                sin(phi1)*sin(phi0))+sin(theta1)*sin(theta0));
            theta_p = tmp*(cos(theta1)*sin(theta0)-sin(theta1)* ...
                cos(theta0)*(cos(phi1)*cos(phi0)+sin(phi1)*sin(phi0)));
            phi_p = tmp*cos(theta1)*cos(theta0)*(cos(phi1)*sin(phi0) ...
                -sin(phi1)*cos(phi0));
            p = [theta_p; phi_p; u1-u0];
            q = [1; 1; 0]; % derivative left (bulk terms apply, see PDEX4)
        end
    end % function surface_fun

    % initial conditions    
    % ------------------------------------------------------------
    function init = ic_fun(zin)
        init = [pretiltl+zin./d_LC*(pretiltr-pretiltl);totaltwist*zin./d_LC; ...
            volts*zin./d_LC];
        % linear guess for u0 isn't good enough, Laplace's eqn must be satisfied
        init(3) = interp1(z, u0, zin);
    end % function ic_fun

    % Simple steady-state solver to obtain the initial distribution of the 
    % electric potential. It is important that u0 is a physical solution, which
    % is a little problematic when alignment layers are present. Although an 
    % analytic solution is available for u0 with a linear tilt profile, it is 
    % rather convoluted. Also need to treat 1 or 2 alignment layer cases 
    % differently. Resort to numerical solution for now, solve the BVP!
    % see u0sol.mw for analytic solution
    function u0 = getu0(z)
        solinit = bvpinit(z,@ic_u);
        sol = bvp4c(@ode_u,@bc_u,solinit);
        u0 = sol.y(1,:);
        function dfdz = ode_u(z,f)
            theta1 = pretiltl+z*(pretiltr-pretiltl);
            theta2 = (pretiltr-pretiltl);
            ezz = epslow + deleps*sin(theta1).^2;
            du2dz = -2*deleps*sin(theta1).*cos(theta1).*theta2.*f(2)./ezz;
            dfdz = [f(2); du2dz];
        end % function ode_u
        function res = bc_u(ya,yb)
            res = [ ya(1)-0;  yb(1)-volts ];
            if d_alignl~=0
                res(1) = (epslow + deleps*sin(pretiltl)^2)*ya(2) ...
                      - eps_alignl*(ya(1)-0)/d_alignl;
            end        
            if d_alignr~=0
                res(2) = (epslow + deleps*sin(pretiltr)^2)*yb(2) ...
                      - eps_alignr*(volts-yb(1))/d_alignr;
            end
        end % function bc_u
        function init = ic_u(z)
               init = [ z*volts; volts ]./d_LC;
        end % function ic_u
    end % function getu0

end % function lc3kdt

% simple Jones code, normal incidence
% might want to interpolate (interp1) onto a finer mesh to improve results
function intensity =  transmittance(z, theta, phi)
    lambda = 0.55;          % wavelength (same units as z)
    no = 1.5; % LC refractive indices (ignore optical effect of align. layers)
    ne = no+lambda/2/max(z); %~halfwave plate;
    pinangle = 45*pi/180;   % polarizer angle
    poutangle = -45*pi/180; % analyzer angle
    Pin  = rot(-pinangle )*[1,0;0,0]*rot(pinangle );
	Pout = rot(-poutangle)*[1,0;0,0]*rot(poutangle);
    d = conv(diff(z), [0.5,0.5]); % thickness of each layer
    J = Pout;
    for iz = 1:length(z),
        J = J*jones(phi(iz), theta(iz), no, ne, d(iz), lambda);
    end
    J = J*Pin;
    intensity = sum(abs(J(:)).^2);
end % function transmittance

function J = jones(twist,tilt,no,ne,d,lambda)
    n = 1/sqrt(sin(tilt)^2/(no*no)+cos(tilt)^2/(ne*ne));
    J = rot(-twist)*[exp(i*2*pi*n*d/lambda),0;0,exp(i*2*pi*no*d/lambda)]*rot(twist);
end % function jones

function R = rot(angle) % rotation matrix
	R = [cos(angle), sin(angle); -sin(angle), cos(angle)];
end % function rot