function [z,theta,phi,u,intensity] = lc3k(voltsin)
%LC3K Find equilibrium liquid crystal orientation in a 1D cell.
%   [Z,THETA,PHI,U,INTENSITY] = LC3K(VOLTS) returns the 
%   tilt (THETA), twist (PHI) and potential (U) as a 
%   function of Z for a 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). Might be advantagous to scale u to lie in
%   [0,1] for the sake for the condition number of the 
%   Jacobian.
%
%   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.
%
%   Another problem is that the initial conditions (set in
%   IC_FUN) must be sufficiently close to the proper
%   solution to obtain the correct result. This function 
%   really needs to be improved to have a good guesses for 
%   Planar, Pi, VAN and VAN cells. Currently, it is set up 
%   to assume a sinusoidal intial tilt. If NUMVOLTSTEPS is 
%   greater than 1, then the voltage is ramped in a number 
%   of iterations and IC_FUN assumes a linear initial tilt.
%   Increase NUMVOLTSTEPS if you do not have a good initial 
%   guess of your own, or if the routine does not converge.
%
%   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 BVP4C, SHOCKBVP.

    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 = 5*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.1;   % left alignment layer thickness  [microns] can be 0
    d_alignr = 0.0;   % right alignment layer thickness [microns] can be 0
    d_LC = 1;%-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 = Inf;%1000e-6; % Weak = ? [J/m^2] or Strong = Inf
    Wr = Inf;%1000e-6; % Weak = ? [J/m^2] or Strong = Inf
    
    % these optimisations only double the speed (since J is complicated), hardly worth the trouble
    options = bvpset('Vectorized', 'on', 'Nmax', 100, 'FJacobian', @jac_3k);

    zinit = linspace(0, d_LC, 11); % initial mesh

    % ramp voltage to obtain solution for high volts or non-planar cells
    numvoltsteps = 1;%ceil(sqrt(e0*(voltsin/d_LC)^2*deleps/K33));
	for voltstep = 1:numvoltsteps,

        volts = voltsin*voltstep^1/numvoltsteps^1;
        disp(sprintf('voltstep %d / %d, %g V', voltstep, numvoltsteps, volts));
        
        solinit = bvpinit(zinit,@ic_fun);
        sol = bvp4c(@ode_3k,@bc_fun,solinit,options);

        z = sol.x; theta = sol.y(1,:); phi = sol.y(3,:); u = sol.y(5,:);%*voltsin;
        intensity = transmittance(z, theta, phi);

        % plot results
        clf reset
        [ax,h1,h2] = plotyy(z,[theta;phi]*180/pi,z,u);
        title(sprintf('Transmittance %g%%', intensity*100));
        set(h1, 'Marker', 'o'); set(h2, 'Marker', 'o');
        xlabel('z')
        ylabel(ax(1), 'tilt, twist [degrees]')
        ylabel(ax(2), ' voltage [V]')
        grid on
        drawnow
	end
    
    % ------------------------------------------------------------
    % minimise Oseen Frank free energy and solve Laplace (vectorized)
    % ------------------------------------------------------------
	function dfdz = ode_3k(z,f)
	% [f(1..6)] = [theta; theta'; phi; phi'; u; u']
        costheta1 = cos(f(1,:));
        sintheta1 = sin(f(1,:));
        theta2 = f(2,:); % the spatial derivative theta'
        phi2 = f(4,:); % the spatial derivative phi'
        u2 = f(6,:); % the spatial derivative u'

        % maple generated (see lc3k.mws)
        dtheta2dz = costheta1.*sintheta1.*(K11*theta2.^2-K33*theta2.^2+ ...
            2*K33*costheta1.^2.*phi2.^2+2*K22*phi2*q0-K33*phi2.^2- ...
            2*K22*costheta1.^2.*phi2.^2-e0*deleps*(u2.*1).^2)./(K11*costheta1.^2 ...
            -K33*costheta1.^2+K33);

        dphi2dz = -2*sintheta1.*theta2.*(2*K33*costheta1.^2.*phi2- ...
            2*K22*costheta1.^2.*phi2+K22*q0-K33*phi2)./costheta1 ...
            ./(-K33*costheta1.^2+K22*costheta1.^2+K33);

        ezz = epslow + deleps*sintheta1.^2;
        du2dz = -2*deleps*sintheta1.*costheta1.*theta2.*u2./ezz;

        dfdz = [theta2; dtheta2dz; phi2; dphi2dz; u2; du2dz];
	end

    % Jacobian function, diff(eultheta3, theta1) etc.
    % doesn't actually help convergence and can't be vectorized either! :(
	function jac = jac_3k(z,f)
        costheta1 = cos(f(1,:));
        sintheta1 = sin(f(1,:));
        theta2 = f(2,:); % the spatial derivative theta'
        phi2 = f(4,:); % the spatial derivative phi'
        u2 = f(6,:); % the spatial derivative u'

        theta2_theta1 = sintheta1^2*(theta2^2*K33-K11*theta2^2-2*K33*costheta1^2 ...
            *phi2^2+2*K22*costheta1^2*phi2^2+e0*deleps*u2^2-2*K22*phi2*q0+K33*phi2^2) ...
            /(-K33*costheta1^2+K11*costheta1^2+K33)-costheta1^2 ...
            *(theta2^2*K33-K11*theta2^2-2*K33*costheta1^2*phi2^2 ...
            +2*K22*costheta1^2*phi2^2+e0*deleps*u2^2-2*K22*phi2*q0+K33 ...
            *phi2^2)/(-K33*costheta1^2+K11*costheta1^2+K33) ...
            -costheta1*sintheta1*(4*K33*costheta1*phi2^2*sintheta1 ...
            -4*K22*costheta1*phi2^2*sintheta1)/(-K33*costheta1^2+K11 ...
            *costheta1^2+K33)+costheta1*sintheta1*(theta2^2*K33-K11* ...
            theta2^2-2*K33*costheta1^2*phi2^2+2*K22*costheta1^2*phi2^2 ...
            +e0*deleps*u2^2-2*K22*phi2*q0+K33*phi2^2)/(-K33*costheta1^2+ ...
            K11*costheta1^2+K33)^2*(2*K33*costheta1*sintheta1 ...
            -2*K11*costheta1*sintheta1);
        theta2_theta2 = -costheta1*sintheta1*(2*theta2*K33-2*K11*theta2)...
            /(-K33*costheta1^2+K11*costheta1^2+K33);
        theta2_phi2 = -costheta1*sintheta1*(-4*K33*costheta1^2*phi2 ...
            +4*K22*costheta1^2*phi2-2*K22*q0+2*K33*phi2)/(-K33*costheta1^2 ...
            +K11*costheta1^2+K33);
        theta2_u2 = -2*costheta1*sintheta1*e0*deleps*u2/(-K33*costheta1^2 ...
            +K11*costheta1^2+K33);

        phi2_theta1 = -2*theta2*(2*K33*costheta1^2*phi2-2*K22*costheta1^2 ...
            *phi2+K22*q0-K33*phi2)/(-K33*costheta1^2+K22*costheta1^2+K33)-2* ...
            sintheta1*theta2*(-4*K33*costheta1*phi2*sintheta1+4*K22 ...
            *costheta1*phi2*sintheta1)/costheta1/(-K33*costheta1^2 ...
            +K22*costheta1^2+K33)-2*sintheta1^2*theta2*(2*K33*costheta1^2 ...
            *phi2-2*K22*costheta1^2*phi2+K22*q0-K33*phi2)/costheta1^2/(-K33* ...
            costheta1^2+K22*costheta1^2+K33)+2*sintheta1*theta2*(2*K33 ...
            *costheta1^2*phi2-2*K22*costheta1^2*phi2+K22*q0-K33*phi2) ...
            /costheta1/(-K33*costheta1^2+K22*costheta1^2+K33)^2*(2*K33 ...
            *costheta1*sintheta1-2*K22*costheta1*sintheta1);
        phi2_theta2 = -2*sintheta1*(2*K33*costheta1^2*phi2-2*K22*costheta1^2 ...
            *phi2+K22*q0-K33*phi2)/costheta1/(-K33*costheta1^2+K22 ...
            *costheta1^2+K33);
        phi2_phi2 = -2*sintheta1*theta2*(2*K33*costheta1^2-2*K22*costheta1^2 ...
            -K33)/costheta1/(-K33*costheta1^2+K22*costheta1^2+K33);

        u2_theta1 = -2*deleps*costheta1^2*theta2*u2/(epslow+deleps*sintheta1^2) ...
            +2*deleps*sintheta1^2*theta2*u2/(epslow+deleps*sintheta1^2) ...
            +4*deleps^2*sintheta1^2*costheta1^2*theta2*u2/(epslow ...
            +deleps*sintheta1^2)^2;
        u2_theta2 = -2*deleps*sintheta1*costheta1*u2/(epslow+deleps*sintheta1^2);
        u2_u2 = -2*deleps*sintheta1*costheta1*theta2/(epslow+deleps*sintheta1^2);

        jac = [	0,1,0,0,0,0; ...
            theta2_theta1, theta2_theta2, 0, theta2_phi2, 0, theta2_u2; ...
            0,0,0,1,0,0; ...
            phi2_theta1, phi2_theta2, 0, phi2_phi2, 0, 0; ...
            0,0,0,0,0,1; ...
            u2_theta1, u2_theta2, 0, 0, 0, u2_u2];
    end % function jac_3k

    % ------------------------------------------------------------
	function res = bc_fun(ya,yb) % boundary conditions

        % note -Wr to take into account the difference in surface normal l&r
        [theta_l, phi_l] = surface_fun(-Wl, ya(1:4), pretiltl, 0); 
        [theta_r, phi_r] = surface_fun( Wr, yb(1:4), pretiltr, totaltwist);
        
        u_l = ya(5)-0; u_r = yb(5)-volts;
        
        % satisfy (z\cdot\varepsilon\cdot u)^+ = (z\cdot\varepsilon\cdot u)^-
        % (i.e. Snell's law) at the interface between the LC and alignment layer       
        if d_alignl~=0
            u_l = (epslow + deleps*sin(pretiltl)^2)*ya(6) ...
                  - eps_alignl*(ya(5)-0)/d_alignl;
        end        
        if d_alignr~=0
            u_r = (epslow + deleps*sin(pretiltr)^2)*yb(6) ...
                  - eps_alignr*(volts-yb(5))/d_alignr;
        end
        
        res = [ theta_l; theta_r; phi_l; phi_r; u_l; u_r ];
    end % function bc_fun

    function [theta_b, phi_b] = surface_fun( W, thetaphi, theta0, phi0) 
    
        theta1 = thetaphi(1); theta2 = thetaphi(2);
        phi1 = thetaphi(3); phi2 = thetaphi(4);
        
        if isinf(W), % strong anchoring, Dirichlet boundary condition
            theta_b = theta1-theta0;
            phi_b = phi1-phi0;
        else % weak anchoring, Neumann boundary condition
            % relate gradients of theta, phi at surface from \deltaF to Fs.
            % boundary terms from the first variation in F
            theta_b = (K11*cos(theta1)^2+K33*sin(theta1)^2)*theta2;
            phi_b = (K22*cos(theta1)^4-K33*cos(theta1)^4+ ...
                K33*cos(theta1)^2)*phi2-K22*q0*cos(theta1)^2;
            % 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_b = theta_b-tmp*(cos(theta1)*sin(theta0)-sin(theta1)* ...
                cos(theta0)*(cos(phi1)*cos(phi0)+sin(phi1)*sin(phi0)));
            phi_b = phi_b-tmp*cos(theta1)*cos(theta0)*( ...
                cos(phi1)*sin(phi0)-sin(phi1)*cos(phi0));
            % anchoring energy is of the form W/2*sin(theta1-theta0)^2
            % theta_b = theta_b+W*1e-6*sin(theta1-theta0)*cos(theta1-theta0);
            % phi_b = phi_b+W*1e-6*sin(phi1-phi0)*cos(phi1-phi0);
        end
    end % function surface_fun
    
    % ------------------------------------------------------------
    function init = ic_fun(z,varargin) % initial conditions for each component
       % of the solution vector, including variables and their derivatives
       if voltstep==1, % linear starting configuration for first iteration
           zs = max(0, min(1, z/d_LC));
           maxtilt = (numvoltsteps~=0)*80*pi/180; % increase for high volts to obtain convergence
           init = [ pretiltl+zs*(pretiltr-pretiltl)+maxtilt*sin(pi*zs);
                    (pretiltr-pretiltl)+maxtilt*pi*cos(pi*zs); ...
                    totaltwist*zs; totaltwist; ...
                    zs*volts; volts];% work with a scaled u in [0,1] % 0.1*volts*zs; 0.1*volts];
       else % use previous solution vector as a initial conditions
           %ind = [1,3:length(sol.x)-2,length(sol.x)]; % assumes alignment layers nodeless
           [dummy, ind] = unique(sol.x);
           init = interp1(sol.x(ind), sol.y(:,(ind))', z)';
       end
    end % function ic_fun

end % function lc3k_Oseen

% 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 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 R = rot(angle) % rotation matrix
	R = [cos(angle), sin(angle); -sin(angle), cos(angle)];
end