#----------------------------------------------------------
# covariance matrix and correlation matrix of sum and difference
# (exercise on the measurements of the sides of the A4 foil)
#
#  GdA 26/2/2016
#------------------------------

m.a = 29.73
m.b = 21.41

sd.a = 0.03    # then change to 0.04 and to 0.05 to see the difference!
sd.b = 0.04
rho.ab = 0.0   # then change to +0.8 to see the difference

cat("\nInput  ---------------------------------------------------------------- \n")
cat(sprintf("a = (%.2f +- %.2f) cm\n", m.a, sd.a ))
cat(sprintf("b = (%.2f +- %.2f) cm\n", m.b, sd.b ))
cat(sprintf("rho(a,b) = %.2f \n", rho.ab ))

#-------------------------------------------------------------------------------
# (unit coefficients made explicit) 

# standard deviations
sd.s = sqrt( sd.a^2 + sd.b^2 + 2*(1)*(1)* (rho.ab*sd.a*sd.b) ) 
sd.d = sqrt(sd.a^2 + sd.b^2 + 2*(1)*(-1)* (rho.ab*sd.a*sd.b) ) 

cat(sprintf("s = (%.2f +- %.2f) cm\n", m.a+m.b, sd.s ))
cat(sprintf("d = (%.2f +- %.2f) cm\n", m.a-m.b, sd.d ))

cov.sd = 1*1*sd.a^2 + 1*(-1)*sd.b^2 +
         ( 1*(-1) + 1*1 ) * ( rho.ab * sd.a * sd.b )  
                                                     
rho.sd = cov.sd / (  sd.s * sd.d )     
                                                                             
cat(sprintf("cov(s,d) = %.3e cm^2\n", cov.sd ))
cat(sprintf("rho(s,d) = %.3f \n", rho.sd))