#Data Electron Microprobe from Nov3, 2005 (1103)
# Use the new paramRot function to calculate the apices of the rotated box
GraphTrue<-c(FALSE, FALSE, FALSE, TRUE)

h1<-w1<-10
h2<-w2<-2

paramRot<- function(typ,xo,yo,xl,yl,th){
	     TH<- th*pi/100
        ROT<-matrix(c(cos(TH),sin(TH),-sin(TH),cos(TH)),2,2)
        XY<-matrix(0,4,2)
	     xl<-xl/2; XY[1:2,1]<-xo-xl; XY[3:4,1]<- xo+xl
	     yl<-yl/2; XY[c(1,4),2]<-yo-yl; XY[2:3,2]<-yo+yl
	     mx<-mean(XY[,1]); my<-mean(XY[,2])
	     MXY<- t(matrix(c(mx,my),2,4))
	     XY<-MXY + t(ROT%*%t(XY-MXY))
	     XY
	     	     }

paramIn<- function(typ,xo,yo,xl,yl,th){
	     TH<- th*pi/100
        ROT<-matrix(c(cos(TH),sin(TH),-sin(TH),cos(TH)),2,2)
        XY<-matrix(0,4,2)
	     xl<-xl/2; XY[1:2,1]<-xo-xl; XY[3:4,1]<- xo+xl
	     yl<-yl/2; XY[c(1,4),2]<-yo-yl; XY[2:3,2]<-yo+yl
	     mx<-mean(XY[,1]); my<-mean(XY[,2])
	     MXY<- t(matrix(c(mx,my),2,4))
	     XY<-MXY + t(ROT%*%t(XY-MXY))
	     
	     yl<-XY[1:2,2]; xl<-XY[1:2,1]
        yr<-XY[3:4,2]; xr<-XY[3:4,1]
        yu<-XY[2:3,2]; xu<-XY[2:3,1]
        yd<-XY[c(4,1),2]; xd<-XY[c(4,1),1]
	     xleft<-lm(xl~yl)
	     xright<-lm(xr~yr)
	     yup<-lm(yu~xu)
	     ydown<-lm(yd~xd)
	     Out<-c(1, xleft[[1]], xright[[1]], ydown[[1]],yup[[1]])
	     Out	
	     }
TestwParam<- function(p,x,y) {
	 In<-FALSE
	 if(p[1]==1) {
	 	         if(p[2]+p[3]*y <= x) {
	 	         	if(p[4]+p[5]*y >= x) {
	 	         		if(p[6]+p[7]*x <= y){
	 	         			if(p[8]+p[9]*x >= y){
	 	         				In<-TRUE}
	 	         			} } }
	 	 	         }
	  In
	  }
	  
Lines<-read.table("Ascii.listarea", header=TRUE)
attach(Lines)
nlin<- length(Lines[,1])
 filcol<-c("red", "green", "blue", "purple", "black", "yellow", "darkgreen")
 fillwd<-c(2    ,  3     ,  2    ,  4      ,  3     ,  3      ,  2    )
filnlev<-c(10    , 6     , 10    ,  6      ,  6     ,  6      ,  10    )
fillty<-c("solid","dashed","solid","solid","dashed","solid", "solid")
filnam<-jpat[1]
# Calcium plots
xy<-read.table(paste(filnam,"_2_Ca.txt",sep=''), skip=33)
attach(xy)
YY<-t(as.matrix(xy))
ny<- length(YY[,1])
nx<- length(YY[1,])
if (GraphTrue[1]==TRUE) {
      quartz(wi=w1, he=h1)
      opar<-par(mar=c(0,0,0,0)+0.05, ann=FALSE)
	   image(YY, col=gray(1:255/255),xlab='', ylab='') 
      mtext(paste(filnam, "Ca++"), line=-1.5, col="white", adj=0.05 )                        
      for (i in 1:nlin) {
       mtext(paste("Xo=", Xo[i], "Yo=",Yo[i], filcol[i]), line=-1.25*(3-((-i)%%3)),  
           col="white", adj= 0.27*ceiling(i/3))
       Pt<-paramRot(1,Xo[i],Yo[i],Xl[i],Yl[i],theta[i])
       xl<-c(Pt[,1],Pt[1,1]); yl<-c(Pt[,2],Pt[1,2])
       lines(xl/nx,yl/ny, col=filcol[i], lwd=2)
                  }  # end for i
                         } # end if (GraphTrue[1]= TRUE)
detach(xy)


# Phosphorous plots
xy<-read.table(paste(filnam,"_3_P.txt",sep=''), skip=33)
attach(xy)
ZZ<-t(as.matrix(xy))
if (GraphTrue[2]== TRUE) {
   quartz(wi=w1, he=w1)
   ppar<-par(mar=c(0,0,0,0)+0.05, ann=FALSE)
   image(ZZ, col=gray(1:255/255),xlab='', ylab='') 
   mtext(paste(filnam, "PO4"), line=-1.5, col="white", adj=0.05 )
   for (i in 1:nlin) {
       mtext(paste("Xo=", Xo[i], "Yo=",Yo[i], filcol[i]), line=-1.25*(3-((-i)%%3)),  
           col="white", adj= 0.27*ceiling(i/3))
    Pt<-paramRot(1,Xo[i],Yo[i],Xl[i],Yl[i],theta[i])
    xl<-c(Pt[,1],Pt[1,1]); yl<-c(Pt[,2],Pt[1,2])
    lines(xl/nx,yl/ny, col=filcol[i], lwd=2)
                  } # end for i #2
                       } # end if (GraphTrue[2]= TRUE)
detach(xy)

MaxCa<- max(YY)
MaxP<- max(ZZ)

# Dot Plots
if (GraphTrue[3]== TRUE) {
 quartz(wi=w1, he=h1)
 qpar<-par(mar=c(0,0,0,0)+0.05, ann=FALSE)
 plot(c(0,1), c(0,1), ty='n')
 for (i in 1:nlin) {
  Pt<-paramIn(1,Xo[i],Yo[i],Xl[i],Yl[i],theta[i])
  for (j in 1:nx) {for (k in 1:ny) {
		TF<-TestwParam(Pt,j,k)
		if(TF==TRUE) {points(ZZ[j,k]/MaxP,YY[j,k]/MaxCa, col=filcol[i], pch=20, cex=0.5)}
		}
	}  # end for j
}  # end for i
	mtext("Ca -->", side = 2, line = -2, cex = 2, adj=0.85)
	mtext("P -->", side = 1, line = -2, cex = 2)
                       } # end if (GraphTrue[3]= TRUE)

# Contour Plots
if (GraphTrue[4]== TRUE) {
  quartz(wi=w1, he=h1)
  rpar<-par(mar=c(0,0,0,0)+0.05, ann=FALSE)
  plot(c(0,1), c(0,1), ty='n')
  for (i in 1:nlin) {
   Pt<-paramIn(1,Xo[i],Yo[i],Xl[i],Yl[i],theta[i])
   dens<-matrix(0,floor(MaxCa/60)+1,floor(MaxP/6)+1)
	
    for (j in 1:nx) {for (k in 1:ny) {
		TF<-TestwParam(Pt,j,k)
		if(TF==TRUE) {dens[floor(YY[j,k]/60)+1,floor(ZZ[j,k]/6)+1]<- dens[floor(YY[j,k]/60)+1,floor(ZZ[j,k]/6)+1] +1}
		}
	}  # end for j
	contour(t(dens)^0.25, add=TRUE, col=filcol[i],
	        nlevels=filnlev[i], lwd=fillwd[i], lty=fillty[i])
}  # end for i
	mtext("Ca -->", side = 2, line = -2, cex = 2, adj=0.85)
	mtext("P -->", side = 1, line = -2, cex = 2)
                           } # end if (GraphTrue[4]= TRUE)


detach(Lines)