%copyright Hongyi Cai, 2007, all rights reserved ;
%contact Hongyi Cai, hongyi.cai@gmail.com ;

%the number of viewing media and the matrix P of all parameters ;
syms k P ;   
k = input('Total number of visual media? \n') ;
P = zeros(k*9, 6) ;


%calculate the on-axis legibility distance D0 ;
syms Ht Rhsw Sd Lb Cp D0 ; 
format short g ;
Sd = input('The denominator of snellen eyesight of the observer? \n') ;
for k = 1:k
    disp('For the visual media of number:') ; 
    disp(k) ;
    Ht = input('The height of text to be viewed in mm? \n') ;
    Rhsw = input('The height-to-strokewidth ratio of the text? \n') ;
    Lb = input('The background luminance level of the text in cd/m2? \n') ;
    Cp = input('The luminance contrast percent of text? \n') ;
    D0 = 2443.5*Ht*(Rhsw^-1)*(Sd^-1)*(Lb^0.213)*(Cp^0.532)/1000 ;  % on-axis viewing distance in meter
    P((9*(k-1)+1):9*k,4) = D0 ; 
end


% The geometries and angles of visual media
syms xc yc zc ;
syms xtl ytl ztl xtm ytm ztm xtr ytr ztr xml yml zml xmr ymr zmr xbl ybl zbl xbm ybm zbm xbr ybr zbr ;
syms h w deltaphi deltaalpha ;

for k = 1:k ;
    disp('For the visual media of number:') ; 
    disp(k) ;
    h = input('height of the visual media in meters? \n') ;
    w = input('width of the visual media in meters? \n') ;
    xc = input('x-coordinate of the center point of visual media in meters? \n') ;
    yc = input('y-coordinate of the center point of visual media in meters? \n') ;
    zc = input('z-coordinate of the center point of visual media in meters? \n') ;
    deltaphi = input('inital horizontal viewing angle of the visual media in degrees? \n') ;
    deltaalpha = input('inital vertical viewing angle of the visual media in degrees? \n') ; 
    
    % other 8 points coordinates in meters calculated from the center point;
    
    % top left point ;
    xtl = xc - (w/2)*cos(deltaphi) + (h/2)*sin(deltaalpha)*sin(deltaphi) ;
    ytl = yc + (w/2)*sin(deltaphi) + (h/2)*sin(deltaalpha)*cos(deltaphi) ;
    ztl = zc + (h/2)*cos(deltaalpha) ;
    
    % top middle point ;
    xtm = xc + (h/2)*sin(deltaalpha)*sin(deltaphi) ;
    ytm = yc + (h/2)*sin(deltaalpha)*cos(deltaphi) ;
    ztm = zc + (h/2)*cos(deltaalpha) ;
    
    % top right point ;
    xtr = xc + (w/2)*cos(deltaphi) + (h/2)*sin(deltaalpha)*sin(deltaphi) ;
    ytr = yc - (w/2)*sin(deltaphi) + (h/2)*sin(deltaalpha)*cos(deltaphi) ;
    ztr = zc + (h/2)*cos(deltaalpha) ;
    
    % middle left point ;
    xml = xc - (w/2)*cos(deltaphi) ;
    yml = yc + (w/2)*sin(deltaphi) ;
    zml = zc ;
    
    % middle right point ;
    xmr = xc + (w/2)*cos(deltaphi) ;
    ymr = yc - (w/2)*sin(deltaphi) ;
    zmr = zc ;
    
    % bottme left point  ;
    xbl = xc - (w/2)*cos(deltaphi) - (h/2)*sin(deltaalpha)*sin(deltaphi) ;
    ybl = yc + (w/2)*sin(deltaphi) - (h/2)*sin(deltaalpha)*cos(deltaphi) ;
    zbl = zc - (h/2)*cos(deltaalpha) ;
    
    % bottom middle point ;
    xbm = xc - (h/2)*sin(deltaalpha)*sin(deltaphi) ;
    ybm = yc - (h/2)*sin(deltaalpha)*cos(deltaphi) ;
    zbm = zc - (h/2)*cos(deltaalpha) ;
    
    % bottom right point ;
    xbr = xc + (w/2)*cos(deltaphi) - (h/2)*sin(deltaalpha)*sin(deltaphi) ;
    ybr = yc - (w/2)*sin(deltaphi) - (h/2)*sin(deltaalpha)*cos(deltaphi) ;
    zbr = zc - (h/2)*cos(deltaalpha) ;
    
    % input all x,y and z coordinates to matrix p ;
    P((9*(k-1) + 1),1) = xtl ;
    P((9*(k-1) + 1),2) = ytl ;     %the coordinates of top left points ;
    P((9*(k-1) + 1),3) = ztl ;
    P((9*(k-1) + 2),1) = xtm ;
    P((9*(k-1) + 2),2) = ytm ;     %the coordinates of top middle points ;
    P((9*(k-1) + 2),3) = ztm ;
    P((9*(k-1) + 3),1) = xtr ;
    P((9*(k-1) + 3),2) = ytr ;     %the coordinates of top right points ;
    P((9*(k-1) + 3),3) = ztr ;
    P((9*(k-1) + 4),1) = xml ;
    P((9*(k-1) + 4),2) = yml ;     %the coordinates of middle left points ;
    P((9*(k-1) + 4),3) = zml ;     
    P((9*(k-1) + 5),1) = xc ;
    P((9*(k-1) + 5),2) = yc ;      %the coordinates of center points ;
    P((9*(k-1) + 5),3) = zc ;
    P((9*(k-1) + 6),1) = xmr ;
    P((9*(k-1) + 6),2) = ymr ;      %the coordinates of middle right points ;
    P((9*(k-1) + 6),3) = zmr ;
    P((9*(k-1) + 7),1) = xbl ;
    P((9*(k-1) + 7),2) = ybl ;      %the coordinates of bottom left points ;
    P((9*(k-1) + 7),3) = zbl ;
    P((9*(k-1) + 8),1) = xbm ;
    P((9*(k-1) + 8),2) = ybm ;      %the coordinates of bottom middle points ;
    P((9*(k-1) + 8),3) = zbm ;
    P((9*(k-1) + 9),1) = xbr ;
    P((9*(k-1) + 9),2) = ybr ;      %the coordinates of bottom right points ;
    P((9*(k-1) + 9),3) = zbr ;
    
    % input deltaphi and deltaalpha into matrix p
    P((9*(k-1)+1):9*k,5) = deltaphi*pi/180 ;  %the delta phi for all 9 points on each visual media ;
    P((9*(k-1)+1):9*k,6) = deltaalpha*pi/180 ; %the delta alpha for all 9 points on each visual media ;    
end


% input the parameters to define the viewing plane
syms seta yv zv ;
seta = input('The angle of sloped viewing plane in degrees? if horizontal viewing plane, input zero. \n') ;
yv = input('The distance in y-coordinate in meters from original point to the start edge of the sloped viewing plane? \n') ;
zv = input('The height in meters of the eyes of the observers on the floor? \n') ;


% start to draw the ideal viewing space of texts presneted on each visual media ;

% define volume -100 meter ~ + 100 meter on x y z coordinates;
[x,y,z] = meshgrid(-100:1:100,-100:1:100,-100:1:100);

% define some angles ;
phi1 = (-90:1:90)*pi/180 ;  %the range of horizontal viewing angle phi ;
alpha1 = (-90:1:90)*pi/180 ;  %the range of vertical viewing angle alpha ;

%define the length of x y z ;
len1 = length(phi1) ;
len2 = length(alpha1) ;
xi = zeros(len1,len2) ;
yi = xi ;
zi = xi ;

% set the values of x, y, z for the contour of center point of first screen ;
for k = 1:k
    for j = 1:9    
        for m = 1:len1
            for n = 1:len2
                phit = phi1(m) ; alphat = alpha1(n) ;  
                    ksi = acos(cos(phit)*cos(alphat))*180/pi ;
                    if (ksi >= 0) && (ksi <= 65.7)
                         xi(m,n) = P((9*(k-1)+j),1) + P((9*(k-1)+j),4)*(cos(phit))^0.5*(cos(alphat))^1.5*sin(phit + P((9*(k-1)+j),5)) ;        
                         yi(m,n) = P((9*(k-1)+j),2) + P((9*(k-1)+j),4)*(cos(phit))^0.5*(cos(alphat))^0.5*(cos(alphat)*cos(phit + P((9*(k-1)+j),5))*cos(P((9*(k-1)+j),6))-sin(alphat)*sin(P((9*(k-1)+j),6)));
                         zi(m,n) = P((9*(k-1)+j),3) + P((9*(k-1)+j),4)*(cos(phit))^0.5*(cos(alphat))^0.5*(cos(alphat)*cos(phit+P((9*(k-1)+j),5))*sin(P((9*(k-1)+j),6))+sin(alphat)*cos(P((9*(k-1)+j),6)));
                     elseif (ksi <= 82.8)
                         xi(m,n) = P((9*(k-1)+j),1) + ((0.024*ksi-0.577)^-1)*P((9*(k-1)+j),4)*(cos(phit))^0.5*(cos(alphat))^1.5*sin(phit + P((9*(k-1)+j),5)) ;        
                         yi(m,n) = P((9*(k-1)+j),2) + ((0.024*ksi-0.577)^-1)*P((9*(k-1)+j),4)*(cos(phit))^0.5*(cos(alphat))^0.5*(cos(alphat)*cos(phit + P((9*(k-1)+j),5))*cos(P((9*(k-1)+j),6))-sin(alphat)*sin(P((9*(k-1)+j),6)));
                         zi(m,n) = P((9*(k-1)+j),3) + ((0.024*ksi-0.577)^-1)*P((9*(k-1)+j),4)*(cos(phit))^0.5*(cos(alphat))^0.5*(cos(alphat)*cos(phit+P((9*(k-1)+j),5))*sin(P((9*(k-1)+j),6))+sin(alphat)*cos(P((9*(k-1)+j),6)));
                     else
                         xi(m,n) = P((9*(k-1)+j),1) ;        
                         yi(m,n) = P((9*(k-1)+j),2) ;
                         zi(m,n) = P((9*(k-1)+j),3) ; 
                    end
              end
         end 
         % plot the contour of center point of the first screen ;
            v = z - y*tan(seta*pi/180)+ yv*tan(seta*pi/180) - zv ;  % Create volume data: defines the plane function
            figure(1) ; hold on ;
            contourslice(x,y,z,v,xi,yi,zi,[0 0]) ;
            view(0,90-seta);
            axis equal ;
            axis ([-50 50 -20 80]) ;
            xlabel('x') ;
            ylabel('y') ;
            zlabel('z') ;
    end
end

% display the parameters ;

disp(' ');
disp ('Results') ;
disp(' ');
disp ('list of the parameters of every 9 viewing points on each visual media');
disp(' ');
disp('          x0           y0           z0          D0          deltaphi      deltaalpha');
disp (P) ;