From 4af0501be6bd5faa170c61c44abc69c653a89527 Mon Sep 17 00:00:00 2001 From: Ilkka Virtanen Date: Fri, 29 Sep 2023 16:35:04 +0300 Subject: [PATCH] Now doing the whole fit in local magnetic coordinates --- R/ISaprioriBAFIM.R | 104 ++++++++++++++++++++++++++++----------------- R/ISfit.3D.R | 11 +++-- R/readPP.3D.R | 54 ++++++++++++++++++++--- 3 files changed, 117 insertions(+), 52 deletions(-) diff --git a/R/ISaprioriBAFIM.R b/R/ISaprioriBAFIM.R index 3442ae9..c07d0cb 100644 --- a/R/ISaprioriBAFIM.R +++ b/R/ISaprioriBAFIM.R @@ -199,6 +199,7 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , PP$param[h,] <- scaleParams(xtmp,PP$apriori[[h]]$parScales,inverse=T) PP$covar[[h]] <- scaleCovar(Stmp,PP$apriori[[h]]$parScales,inverse=T) PP$std[h,] <- sqrt(diag(PP$covar[[h]])) +# print(height[h]) } } @@ -220,19 +221,20 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , # The correlation powers solved from known variances, height steps, and correlation lengths corrP <- PP$std[,1:12]**2 - corrP[,1] <- corrP[,1]*dheights/(BAFIMpar$Ne[3]*hsAlt) # Ne - corrP[,2] <- corrP[,2]*dheights/(BAFIMpar$Ti[3]*hsAlt) # Tipar - corrP[,3] <- corrP[,3]*dheights/(BAFIMpar$Ti[3]*hsAlt) # Tiperp - corrP[,4] <- corrP[,4]*dheights/(BAFIMpar$Te[3]*hsAlt) # Tepar - corrP[,5] <- corrP[,5]*dheights/(BAFIMpar$Te[3]*hsAlt) # Teperp - corrP[,6] <- corrP[,6]*dheights/(BAFIMpar$Coll[3]*hsAlt) # Collisions - corrP[,7] <- corrP[,7]*dheights/(BAFIMpar$Viperp[3]*hsAlt) # Vix - corrP[,8] <- corrP[,8]*dheights/(BAFIMpar$Viperp[3]*hsAlt) # Viy - corrP[,9] <- corrP[,9]*dheights/(BAFIMpar$Vipar[3]*hsAlt) # ViB - corrP[,10] <- corrP[,10]*dheights/(BAFIMpar$Mp[3]*hsAlt) # Molecular ions - corrP[,11] <- corrP[,11]*dheights/(BAFIMpar$Op[3]*hsAlt) # O+ - corrP[,12] <- corrP[,12]*dheights/(BAFIMpar$Hp[3]*hsAlt) # H+ - + corrP[,1] <- corrP[,1] * dheights / (BAFIMpar$Ne[3]*hsAlt) # Ne + corrP[,2] <- corrP[,2] * dheights / (BAFIMpar$Ti[3]*hsAlt) # Tipar + corrP[,3] <- corrP[,3] * dheights / (BAFIMpar$Ti[3]*hsAlt) # Tiperp + corrP[,4] <- corrP[,4] * dheights / (BAFIMpar$Te[3]*hsAlt) # Tepar + corrP[,5] <- corrP[,5] * dheights / (BAFIMpar$Te[3]*hsAlt) # Teperp + corrP[,6] <- corrP[,6] * dheights / (BAFIMpar$Coll[3]*hsAlt) # Collisions + corrP[,7] <- corrP[,7] * dheights / (BAFIMpar$Viperp[3]*hsAlt) # Vix + corrP[,8] <- corrP[,8] * dheights / (BAFIMpar$Viperp[3]*hsAlt) # Viy + corrP[,9] <- corrP[,9] * dheights / (BAFIMpar$Vipar[3]*hsAlt) # ViB + corrP[,10] <- corrP[,10] * dheights / (BAFIMpar$Mp[3]*hsAlt) # Molecular ions + corrP[,11] <- corrP[,11] * dheights / (BAFIMpar$Op[3]*hsAlt) # O+ + corrP[,12] <- corrP[,12] * dheights / (BAFIMpar$Hp[3]*hsAlt) # H+ + + #The first order terms. # M is always zero for the first and higher order terms @@ -264,18 +266,18 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , A[Aind,hind] <- -2 A[Aind,hind+1] <- 1 - SNe[Aind] <- 8 * corrP[hind,1] * ( dheights[hind] / (BAFIMpar$Ne[3] * hsAlt[hind]) )**3 - STipar[Aind] <- 8 * corrP[hind,2] * ( dheights[hind] / (BAFIMpar$Ti[3] * hsAlt[hind]) )**3 - STiperp[Aind] <- 8 * corrP[hind,3] * ( dheights[hind] / (BAFIMpar$Ti[3] * hsAlt[hind]) )**3 - STepar[Aind] <- 8 * corrP[hind,4] * ( dheights[hind] / (BAFIMpar$Te[3] * hsAlt[hind]) )**3 - STeperp[Aind] <- 8 * corrP[hind,5] * ( dheights[hind] / (BAFIMpar$Te[3] * hsAlt[hind]) )**3 - SColl[Aind] <- 8 * corrP[hind,6] * ( dheights[hind] / (BAFIMpar$Coll[3] * hsAlt[hind]) )**3 - SVix[Aind] <- 8 * corrP[hind,7] * ( dheights[hind] / (BAFIMpar$Viperp[3] * hsAlt[hind]) )**3 - SViy[Aind] <- 8 * corrP[hind,8] * ( dheights[hind] / (BAFIMpar$Viperp[3] * hsAlt[hind]) )**3 - SVipar[Aind] <- 8 * corrP[hind,9] * ( dheights[hind] / (BAFIMpar$Vipar[3] * hsAlt[hind]) )**3 - SMp[Aind] <- 8 * corrP[hind,10] * ( dheights[hind] / (BAFIMpar$Mp[3] * hsAlt[hind]) )**3 - SOp[Aind] <- 8 * corrP[hind,11] * ( dheights[hind] / (BAFIMpar$Op[3] * hsAlt[hind]) )**3 - SHp[Aind] <- 8 * corrP[hind,12] * ( dheights[hind] / (BAFIMpar$Hp[3] * hsAlt[hind]) )**3 + SNe[Aind] <- 8 * corrP[hind,1] * ( dheights[hind] / (BAFIMpar$Ne[3] * hsAlt[hind]) )**3 + STipar[Aind] <- 8 * corrP[hind,2] * ( dheights[hind] / (BAFIMpar$Ti[3] * hsAlt[hind]) )**3 + STiperp[Aind] <- 8 * corrP[hind,3] * ( dheights[hind] / (BAFIMpar$Ti[3] * hsAlt[hind]) )**3 + STepar[Aind] <- 8 * corrP[hind,4] * ( dheights[hind] / (BAFIMpar$Te[3] * hsAlt[hind]) )**3 + STeperp[Aind] <- 8 * corrP[hind,5] * ( dheights[hind] / (BAFIMpar$Te[3] * hsAlt[hind]) )**3 + SColl[Aind] <- 8 * corrP[hind,6] * ( dheights[hind] / (BAFIMpar$Coll[3] * hsAlt[hind]) )**3 + SVix[Aind] <- 8 * corrP[hind,7] * ( dheights[hind] / (BAFIMpar$Viperp[3] * hsAlt[hind]) )**3 + SViy[Aind] <- 8 * corrP[hind,8] * ( dheights[hind] / (BAFIMpar$Viperp[3] * hsAlt[hind]) )**3 + SVipar[Aind] <- 8 * corrP[hind,9] * ( dheights[hind] / (BAFIMpar$Vipar[3] * hsAlt[hind]) )**3 + SMp[Aind] <- 8 * corrP[hind,10] * ( dheights[hind] / (BAFIMpar$Mp[3] * hsAlt[hind]) )**3 + SOp[Aind] <- 8 * corrP[hind,11] * ( dheights[hind] / (BAFIMpar$Op[3] * hsAlt[hind]) )**3 + SHp[Aind] <- 8 * corrP[hind,12] * ( dheights[hind] / (BAFIMpar$Hp[3] * hsAlt[hind]) )**3 Aind <- Aind + 1 } @@ -358,18 +360,18 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , HpCorr <- Xpost[(11*nh+1):(12*nh)]; #Standard deviations. NOTICE: we could pick the full covariance matrices at each height! - NeErrCorr <- sqrt(diag(Cpost[1:nh,1:nh])); - TiparErrCorr <- sqrt(diag(Cpost[(nh+1):(2*nh),(nh+1):(2*nh)])); - TiperpErrCorr <- sqrt(diag(Cpost[(2*nh+1):(3*nh),(2*nh+1):(3*nh)])); - TeparErrCorr <- sqrt(diag(Cpost[(3*nh+1):(4*nh),(3*nh+1):(4*nh)])); - TeperpErrCorr <- sqrt(diag(Cpost[(4*nh+1):(5*nh),(4*nh+1):(5*nh)])); - CollErrCorr <- sqrt(diag(Cpost[(5*nh+1):(6*nh),(5*nh+1):(6*nh)])); - VixErrCorr <- sqrt(diag(Cpost[(6*nh+1):(7*nh),(6*nh+1):(7*nh)])); - ViyErrCorr <- sqrt(diag(Cpost[(7*nh+1):(8*nh),(7*nh+1):(8*nh)])); - ViparErrCorr <- sqrt(diag(Cpost[(8*nh+1):(9*nh),(8*nh+1):(9*nh)])); - MpErrCorr <- sqrt(diag(Cpost[(9*nh+1):(10*nh),(9*nh+1):(10*nh)])); - OpErrCorr <- sqrt(diag(Cpost[(10*nh+1):(11*nh),(10*nh+1):(11*nh)])); - HpErrCorr <- sqrt(diag(Cpost[(11*nh+1):(12*nh),(11*nh+1):(12*nh)])); + NeErrCorr <- sqrt(diag(Cpost[ 1:nh, 1:nh])); + TiparErrCorr <- sqrt(diag(Cpost[( nh+1): (2*nh), (nh+1): (2*nh)])); + TiperpErrCorr <- sqrt(diag(Cpost[( 2*nh+1): (3*nh), (2*nh+1): (3*nh)])); + TeparErrCorr <- sqrt(diag(Cpost[( 3*nh+1): (4*nh), (3*nh+1): (4*nh)])); + TeperpErrCorr <- sqrt(diag(Cpost[( 4*nh+1): (5*nh), (4*nh+1): (5*nh)])); + CollErrCorr <- sqrt(diag(Cpost[( 5*nh+1): (6*nh), (5*nh+1): (6*nh)])); + VixErrCorr <- sqrt(diag(Cpost[( 6*nh+1): (7*nh), (6*nh+1): (7*nh)])); + ViyErrCorr <- sqrt(diag(Cpost[( 7*nh+1): (8*nh), (7*nh+1): (8*nh)])); + ViparErrCorr <- sqrt(diag(Cpost[( 8*nh+1): (9*nh), (8*nh+1): (9*nh)])); + MpErrCorr <- sqrt(diag(Cpost[( 9*nh+1) :(10*nh), (9*nh+1):(10*nh)])); + OpErrCorr <- sqrt(diag(Cpost[(10*nh+1) :(11*nh),(10*nh+1):(11*nh)])); + HpErrCorr <- sqrt(diag(Cpost[(11*nh+1) :(12*nh),(11*nh+1):(12*nh)])); # plot(Xpost[(3*nh+1):(4*nh)],height,xlim=c(0,3000)) # lines(PP$param[,4],height) @@ -389,63 +391,75 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , nesmooth[nesmooth<=1] <- 1 lines((nesmooth),height) lines((nesmooth+NeErrCorr),height,col='red') + lines((nesmooth+NeErrCorr+BAFIMpar$Ne[4]*sqrt(dt)),height,col='green') plot(PP$param[,2],height,xlim=c(0,3000)) lines(TiparCorr,height) lines(PP$param[,2]+PP$std[,2],height,col='blue') lines(TiparCorr+TiparErrCorr,height,col='red') + lines((TiparCorr+TiparErrCorr+BAFIMpar$Ti[4]*sqrt(dt)),height,col='green') plot(PP$param[,3],height,xlim=c(0,3000)) lines(TiperpCorr,height) lines(PP$param[,3]+PP$std[,3],height,col='blue') lines(TiperpCorr+TiperpErrCorr,height,col='red') + lines((TiperpCorr+TiperpErrCorr+BAFIMpar$Ti[4]*sqrt(dt)),height,col='green') plot(PP$param[,4],height,xlim=c(0,3000)) lines(TeparCorr,height) lines(PP$param[,4]+PP$std[,4],height,col='blue') lines(TeparCorr+TeparErrCorr,height,col='red') + lines((TeparCorr+TeparErrCorr+BAFIMpar$Te[4]*sqrt(dt)),height,col='green') plot(PP$param[,5],height,xlim=c(0,3000)) lines(TeperpCorr,height) lines(PP$param[,5]+PP$std[,5],height,col='blue') lines(TeperpCorr+TeperpErrCorr,height,col='red') + lines((TeperpCorr+TeperpErrCorr+BAFIMpar$Te[4]*sqrt(dt)),height,col='green') plot(PP$param[,6],height,xlim=c(0,1e5)) lines(CollCorr,height) lines(PP$param[,6]+PP$std[,6],height,col='blue') lines(CollCorr+CollErrCorr,height,col='red') + lines((CollCorr+CollErrCorr+BAFIMpar$Coll[4]*sqrt(dt)),height,col='green') plot(PP$param[,7],height,xlim=c(-1,1)*100) lines(VixCorr,height) lines(PP$param[,7]+PP$std[,7],height,col='blue') lines(VixCorr+VixErrCorr,height,col='red') + lines((VixCorr+VixErrCorr+BAFIMpar$Viperp[4]*sqrt(dt)),height,col='green') plot(PP$param[,8],height,xlim=c(-1,1)*100) lines(ViyCorr,height) lines(PP$param[,8]+PP$std[,8],height,col='blue') lines(ViyCorr+ViyErrCorr,height,col='red') + lines((ViyCorr+ViyErrCorr+BAFIMpar$Viperp[4]*sqrt(dt)),height,col='green') plot(PP$param[,9],height,xlim=c(-1,1)*100) lines(ViparCorr,height) lines(PP$param[,9]+PP$std[,9],height,col='blue') lines(ViparCorr+ViparErrCorr,height,col='red') + lines((ViparCorr+ViparErrCorr+BAFIMpar$Vipar[4]*sqrt(dt)),height,col='green') plot(PP$param[,10],height,xlim=c(0,1)) lines(MpCorr,height) lines(PP$param[,10]+PP$std[,10],height,col='blue') lines(MpCorr+MpErrCorr,height,col='red') + lines((MpCorr+MpErrCorr+BAFIMpar$Mp[4]*sqrt(dt)),height,col='green') plot(PP$param[,11],height,xlim=c(0,1)) lines(OpCorr,height) lines(PP$param[,11]+PP$std[,11],height,col='blue') lines(OpCorr+OpErrCorr,height,col='red') + lines((OpCorr+OpErrCorr+BAFIMpar$Op[4]*sqrt(dt)),height,col='green') plot(PP$param[,12],height,xlim=c(0,1)) lines(HpCorr,height) lines(PP$param[,12]+PP$std[,12],height,col='blue') lines(HpCorr+HpErrCorr,height,col='red') + lines((HpCorr+HpErrCorr+BAFIMpar$Hp[4]*sqrt(dt)),height,col='green') # dt <- abs(as.double(ISOdate(date[1],date[2],date[3],date[4],date[5],date[6])) - PP$time_sec) @@ -508,7 +522,7 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , nPar <- length(aprioriParam) # number of imaginary apriori "measurements" - nApriori <- ( nPar + 5 ) + nApriori <- ( nPar + 6 ) # apriori theory matrix aprioriTheory <- matrix( 0 , nrow=nApriori , ncol=nPar ) @@ -803,6 +817,11 @@ ISaprioriBAFIM <- function( PP , date , latitude , longitude , height , nSite , aprioriMeas[c(4,5)] <- 0 } curRow <- curRow + 1 + aprioriTheory[curRow,c(3,5)] <- c(1,-1) + aprioriMeas[curRow] <- 0 + aprioriStd[curRow] <- ifelse(height[h]nburnin,TRUE,FALSE) , ... ) @@ -536,9 +538,6 @@ ISfit.3D <- function( ddirs='.' , odir='.' , heightLimits.km=NA , timeRes.s=60 # scale the results back to physical units and add some metadata for(h in seq(nh)){ if(fitGate[h]){ - # rotate the velocities back to geodetic coordinate system (for backward compatibility) - BrotInv <- solve(Brotmat[[h]]) - param[h,7:9] <- param[h,7:9]%*%BrotInv param[h,] <- scaleParams( fitpar[[h]]$param , scale=apriori[[h]]$parScales , inverse=TRUE ) covar[[h]] <- scaleCovar( fitpar[[h]]$covar , scale=apriori[[h]]$parScales , inverse=TRUE) std[h,] <- sqrt(diag(covar[[h]])) @@ -576,7 +575,7 @@ ISfit.3D <- function( ddirs='.' , odir='.' , heightLimits.km=NA , timeRes.s=60 # save the results to file - PP <- list(param=param,std=std,model=model,chisqr=chisqr,status=status,time_sec=time_sec,date=date,POSIXtime=POSIXtime,height=height,latitude=latitude,longitude=longitude,sites=sites,intersect=intersect,covar=covar,B=B,heightLimits.km=hlims/1000,contribSites=contribSites,mIon=c(30.5,16.0,1.0),MCMC=MCMC,timeLimits.s=iperLimits[k:(k+1)],functionCall=functionCall,apriori=apriori,resFile=resFile , resDir=odir) + PP <- list(param=param,std=std,model=model,chisqr=chisqr,status=status,time_sec=time_sec,date=date,POSIXtime=POSIXtime,height=height,latitude=latitude,longitude=longitude,sites=sites,intersect=intersect,covar=covar,B=B,heightLimits.km=hlims/1000,contribSites=contribSites,mIon=c(30.5,16.0,1.0),MCMC=MCMC,timeLimits.s=iperLimits[k:(k+1)],functionCall=functionCall,apriori=apriori,resFile=resFile , resDir=odir,ViCoordinates='ENUmagnetic') if(nnn>nburnin){ save( PP , file=file.path(odir,resFile) ) } diff --git a/R/readPP.3D.R b/R/readPP.3D.R index 1f67bb1..0d70a2f 100644 --- a/R/readPP.3D.R +++ b/R/readPP.3D.R @@ -41,6 +41,11 @@ readPP.3D <- function(dpath,measuredOnly=T,nSiteVi=3,recursive=F,...){ # number of receiver sites nSites <- dim(PP$sites)[1] + # coordinate system for Vi (used to be geographic ENU, but geomagnetic in later versions) + if (is.null( ViCoord <- PP$ViCoordinates ) ){ + ViCoor <- 'ENUgeodetic' + } + # allocate the necessary arrays param <- array(NA,dim=c(nHeight,nPar+4*nSites+7,nFile)) std <- array(NA,dim=c(nHeight,nPar+4*nSites+7,nFile)) @@ -136,8 +141,16 @@ readPP.3D <- function(dpath,measuredOnly=T,nSiteVi=3,recursive=F,...){ param[r,nPar+1,k] <- (PP$param[r,2] + 2*PP$param[r,3] ) / 3 # electron temperature (Te_par + 2*Te_perp)/3 param[r,nPar+2,k] <- (PP$param[r,4] + 2*PP$param[r,5] ) / 3 + # ion velocity along magnetic field (positive upward) - param[r,nPar+3,k] <- -PP$param[r,7:9]%*%PP$B[r,]/sqrt(sum(PP$B[r,]**2)) + if(ViCoord=='ENUgeodetic'){ + param[r,nPar+3,k] <- -PP$param[r,7:9]%*%PP$B[r,]/sqrt(sum(PP$B[r,]**2)) + }else if(ViCoord=='ENUgeomagnetic'){ + param[r,nPar+3,k] <- PP$param[r,7:9] + }else{ + stop('Unknown coordinate system for Vi') + } + # ion perpendicular/parallel temperature ratio # approximation of the ratio of two correlated normal random variables # from Hayya et al., A note on the ratio of two normally distributed variables, @@ -162,7 +175,13 @@ readPP.3D <- function(dpath,measuredOnly=T,nSiteVi=3,recursive=F,...){ std[r,nPar+1,k] <- sqrt(c(1,2)%*%PP$covar[[r]][2:3,2:3]%*%c(1,2)/9) std[r,nPar+2,k] <- sqrt(c(1,2)%*%PP$covar[[r]][4:5,4:5]%*%c(1,2)/9) - std[r,nPar+3,k] <- sqrt(PP$B[r,]%*%PP$covar[[r]][7:9,7:9]%*%PP$B[r,]/sum(PP$B[r,]**2)) + if(ViCoord=='ENUgeodetic'){ + std[r,nPar+3,k] <- sqrt(PP$B[r,]%*%PP$covar[[r]][7:9,7:9]%*%PP$B[r,]/sum(PP$B[r,]**2)) + }else if(ViCoord=='ENUgeomagnetic'){ + std[r,nPar+3,k] <- sqrt(diag(PP$covar[[r]][7:9,7:9])) + }else{ + stop('Unknown coordinate system for Vi') + } # std[r,nPar+4,k] <- sqrt( PP$covar[[r]][2,2]*PP$param[r,3]**2/PP$param[r,2]**4 + # PP$covar[[r]][3,3]/PP$param[r,2]**2 #- ## 2*PP$covar[[r]][2,3]*PP$param[r,3]/PP$param[r,2]**3 @@ -173,7 +192,7 @@ readPP.3D <- function(dpath,measuredOnly=T,nSiteVi=3,recursive=F,...){ # 2*PP$covar[[r]][4,5]*PP$param[r,5]/PP$param[r,4]**3 # ) std[r,nPar+5,k] <- sqrt(PP$covar[[r]][4,4]+PP$covar[[r]][5,5]+2*PP$covar[[r]][4,5]) - # ion velocity components perpendicular to the magnetic field.. + # ion velocity components perpendicular to the magnetic field, or in ENU if ViCooord==ENUgeomagnetic # geomagnetic north Bhor <- c(PP$B[r,1:2],0) # geomagnetic east is perpendicular both to B and Bhor @@ -183,10 +202,31 @@ readPP.3D <- function(dpath,measuredOnly=T,nSiteVi=3,recursive=F,...){ By <- radarPointings:::vectorProduct.cartesian(Bx,PP$B[r,]) By <- By / sqrt(sum(By**2)) - param[r,nPar+6,k] <- PP$param[r,7:9]%*%Bx - param[r,nPar+7,k] <- PP$param[r,7:9]%*%By - std[r,nPar+6,k] <- sqrt(Bx%*%PP$covar[[r]][7:9,7:9]%*%Bx) - std[r,nPar+7,k] <- sqrt(By%*%PP$covar[[r]][7:9,7:9]%*%By) + if(ViCoord=='ENUgeodetic'){ + param[r,nPar+6,k] <- PP$param[r,7:9]%*%Bx + param[r,nPar+7,k] <- PP$param[r,7:9]%*%By + std[r,nPar+6,k] <- sqrt(Bx%*%PP$covar[[r]][7:9,7:9]%*%Bx) + std[r,nPar+7,k] <- sqrt(By%*%PP$covar[[r]][7:9,7:9]%*%By) + }else if(ViCoord=='ENUgeomagnetic'){ + + rotmat <- matrix(c(Bx,By,Bz),byrow=F,ncol=3) + irotmat <- solve(rotmat) + + # put the Vi in geomagnetic and geodetic systems in the same places as above to avoid confusion + param[r,nPar+6,k] <- PP$param[r,7] + param[r,7,k] <- PP$param[r,7:9]%*%irotmat[,1] + param[r,nPar+7,k] <- PP$param[r,8] + param[r,8,k] <- PP$param[r,7:9]%*%irotmat[,2] + param[r,9,k] <- PP$param[r,7:9]%*%irotmat[,3] + + std[r,nPar+6,k] <- sqrt(PP$covar[[r]][7,7]) + std[r,7,k] <- sqrt(irotmat[,1]%*%PP$covar[[r]][7:9,7:9]%*%irotmat[,1]) + std[r,nPar+7,k] <- sqrt(PP$covar[[r]][8,8]) + std[r,8,k] <- sqrt(irotmat[,2]%*%PP$covar[[r]][7:9,7:9]%*%irotmat[,2]) + std[r,9,k] <- sqrt(irotmat[,3]%*%PP$covar[[r]][7:9,7:9]%*%irotmat[,3]) + }else{ + stop('Unknown coordinate system for Vi') + } for( s in PP$contribSites[[r]]){