prog='JPDF'
%----------------------------------------------------------------------------------------
%          Visualisation of Interaction between two independent Variables
%          ==============================================================
%
% DESCRIPTION:
% ------------
% This program calculates and visualises the joint probability between two
% independent variable. Since the here treated topics fall in the area of
% Wind Engineering on variable will quite likely be the wind. In priciple
% both variables can individually be defined.
%
% Program ID:
% -----------
% File name       : JPDF.m
% Author          : hko
% Development Log : 2011-01-24  hko   Basic structure of the program
%                   2013-03-14  hko   Visualisation of JPDF in four graphs
%                                     Decision function and JP over decision area.
%
%
%========================================================================================
 close all
 clear all
%========================================================================================
% 1) DEFINITIONS AND CENTRAL ADJUSTMENTS OF THE PROGRAM:
%    ---------------------------------------------------
% 1.1) Description of Text for Plots:
%      ------------------------------------------

pi     = 4*atan(1.);

DoPDF  = 1; % switch for plotting the PROJECTED shape of PDFs for both
            % variables on the side walls of the 3D graph.

% 1.2) Distribution densities of variables:
%      ------------------------------------
%      M1, M2 = mean value of both variables
%      S1, S2 = standard deviation of both variables

% Variable 1: Daily Mean Wind Speeds (10-minutes mean)
%             2-parametric Weibull Distribution

A = 5;
k = 2;

% Variable 2: Air Temperatures (degC)
%             Normal Distribution

M2 = 15;
S2 = 6;


%========================================================================================
% 2) CONTRUCTION OF DISTRIBUTION DENSITIES:
%    --------------------------------------
% 2.1 Definition of calculation settings

NU = 100;   % discretisation of the velocity axis 0-30m/s in 0.3m/s steps
NT = 100;   % discretisation of temperature axis -10 to 40degC in 0.5degC steps
dU = 0.2;   % wind velocity step width [m/s]
dT = 0.5;   % air temperature step width [degC]

U0 = 0;     % lowest wind speed [m/s]
T0 = -10;   % lowest air temperature [degC]

U    = zeros(1,NU);  % vector for wind speed range
T    = zeros(1,NT);  % vector for airtemperature range
R    = zeros(NT,NU); % Result matrix
D    = zeros(NT,NU); % Decision matrix
pdfU = zeros(1,NU);  % PDF vector for wind velocity
pdfT = zeros(1,NT);  % PDF vector for air temperature
cdfU = zeros(1,NU);  % CDF vector for wind velocity
cdfT = zeros(1,NT);  % CDF vector for air temperature

% 2.2 Probility Densities Functions (PDF)of individual Variables
% --------------------------------------------------------------

for i=1:NU   % Wind Density Distribution
    U(i)    = U0+(i-1)*dU;
    pdfU(i) = k/A*(U(i)/A)^(k-1)*exp(-(U(i)/A)^k);
end

cdfU(1) = pdfU(1)*dU;
for i=2:NU
    cdfU(i)=pdfU(i)*dU+cdfU(i-1);
end

for i=1:NT   % Temperature Density Distribution
    T(i)    = T0+(i-1)*dT;
    pdfT(i) = 1/(S2*sqrt(2*pi))*exp(-0.5*((T(i)-M2)/S2)^2);
end

cdfT(1) = pdfT(1)*dT;
for i=2:NT
    cdfT(i)=pdfT(i)*dT+cdfT(i-1);
end

% 2.3 Joint Probability Matrix R
% ------------------------------

for i=1:NU
   for j=1:NT
      R(j,i)=pdfU(i)*pdfT(j);
   end
end

% 2.4 Decision function and probability of decision area
% ------------------------------------------------------

JP  = 0; % Value of jount probability in decision area
JPt = 0; % Value of total joint probability

for i=1:NU
   for j=1:NT
      JPt = JPt+R(j,i)*dU*dT;
      if (U(i)>=6)&&(T(j)<=2)       % DECISION FUNCTION
         D(j,i) = 1.;
         JP     = JP+R(j,i)*dU*dT;
      end
   end
end

% Check for probability beyond decision point in individual PDFs:

pU6=0;
for i=1:NU
    if U(i)>=6;pU6=pU6+pdfU(i)*dU;end
end

pT2=0;
for i=1:NT
    if T(i)<=2;pT2=pT2+pdfT(i)*dT;end
end

fprintf(1,'Probability of u>=6m/s                   : %7.5f [-]\n',pU6)
fprintf(1,'Probability of T<=2degC                  : %7.5f [-]\n',pT2)
fprintf(1,'Joint probability                        : %7.5f [-]\n',pU6*pT2)
fprintf(1,' \n')
fprintf(1,'Value of total joint probability         : %7.5f [-]\n',JPt)
fprintf(1,'Value of decision space joint probability: %7.5f [-]\n',JP)

        
%========================================================================================
% 3) GRAPHICAL DISPLAY OF THE EXTREME VALUE ANALYSIS:
%----------------------------------------------------
% Display Definitions
% -------------------

scrsz = get(0,'ScreenSize');

%========================================================================================

figure('Name','3D Joint Probability Density','Position',[5 0.40*scrsz(4) 0.5*scrsz(3) 0.5*scrsz(4)])

R(1,1)=-0.000001; % This point in the joint probability matrix is assigned 
                  % artificially with a negative value to activate the 
                  % offset for the isolines below the 3D graph.
surfc(U,T,R,'EdgeColor','none');hold on

if DoPDF==1
    X1 = zeros(1,NT); X1 = X1+U(NU);
    Y1 = zeros(1,NU); Y1 = Y1+T(NT);
    fact1 = max(pdfU);
    fact2 = max(pdfT);
    plot3(X1,T,pdfT*fact1,'--r','MarkerSize',1); hold on
    plot3(U,Y1,pdfU*fact2,'--b','MarkerSize',1); hold on
end


view(-38,18);
title('Joint Probability Density');
xlabel('wind speed [m/s]');
ylabel('air temperature [degC]');
zlabel('probability of occurrence [-]');

eval(['print -dtiff -zbuffer JointProbability3D']);

%----------------------------------------------------------------------------------------

figure('Name','Probability Isolines','Position',[0.515*scrsz(3) 0.4*scrsz(4) 0.48*scrsz(3) 0.5*scrsz(4)])
contour(U,T,R,20)
hold on
for i=1:NU
   for j=1:NT
      if D(j,i)==1
         plot(U(i),T(j),'.b','MarkerSize',2); hold on
      end
   end
end

title('Isolines of Joint Probability');
xlabel('wind speed [m/s]');
ylabel('air temperature [degC]');
grid on

eval(['print -dtiff -zbuffer JPDisolines']);

%----------------------------------------------------------------------------------------

figure('Name','Distribution Density: Wind Speed','Position',[5 0.06*scrsz(4) 0.3*scrsz(3) 0.24*scrsz(4)])
plot(U,pdfU,'-b','LineWidth',2)
xlabel('wind speed [m/s]');
ylabel('probability of occurrence [-]');
grid on

eval(['print -dtiff -zbuffer pdfX1']);

%----------------------------------------------------------------------------------------

figure('Name','Distribution Density: Air temperature','Position',[0.316*scrsz(3) 0.06*scrsz(4) 0.3*scrsz(3) 0.24*scrsz(4)])
plot(T,pdfT,'-r','LineWidth',2)
xlabel('air temperature [degC]');
ylabel('probability of occurrence [-]');
grid on

eval(['print -dtiff -zbuffer pdfX2']);