ROTGEN2<-function(plat,plong,glat,glong) {
#generates 3 by 3 matrix - rotates standard coordinate system to new one
#pole (plat,plong) is rotated to North pole = (0,0,1)
#"Greenwich" (glat,glong) is rotated to a point with longitude 0
eta<-.Fortran("trans1",plat,plong,rep(1.2,3))[[3]]
xcenter<-.Fortran("trans1",glat,glong,rep(1.2,3))[[3]]
xcenter<-xcenter - sum(xcenter*eta)*eta
xcenter<-xcenter/sqrt(sum(xcenter^2))
x2<-c(eta[2]*xcenter[3]-eta[3]*xcenter[2],
      eta[3]*xcenter[1]-eta[1]*xcenter[3],
      eta[1]*xcenter[2]-eta[2]*xcenter[1])
etamat<-cbind(xcenter,x2,eta)
t(etamat)
}

CONVMAT2<-function(n, xmat) {
#function to change a matrix stored in full mode to symmetric mode
        xsim <- NULL
        for(i in 1:n) {
                for(j in 1:i)
                        xsim <- c(xsim, xmat[i, j])
        }
}

BINGHAM<-function(fitrot,df1=3,alpha=.95,indside=c(1,2),plat=90,plong=0,glat=0,glong=0) {
#calls bingham subroutine
#rotation is from plate indside[1] to indside[2]
call<-sys.call()
qhat<-as.double(1:4)
if (length(fitrot$names)>0) pnames<-fitrot$names[indside] else pnames<-NULL

if (length(fitrot$alat)==3) {
#1 ridge from a triple junction
  xmat<-fitrot$Xmat
  xtx<-t(xmat)%*%xmat
  invxtx<-solve(xtx)
  if ((indside[1]==1)&(indside[2]==2)) {alat<-fitrot$alat[1]
    along<-fitrot$along[1]
    rho<-fitrot$rho[1]
    H11.2<-solve(invxtx[1:3,1:3])
    ahat<-fitrot$ahat12}
  if ((indside[1]==2)&(indside[2]==1)) {alat<- -fitrot$alat[1]
    along<-fitrot$along[1]+180
    rho<-fitrot$rho[1]
    H11.2<-fitrot$ahat12%*%solve(invxtx[1:3,1:3])%*%t(fitrot$ahat12)
    ahat<-t(fitrot$ahat12)}
  if ((indside[1]==1)&(indside[2]==3)) {alat<-fitrot$alat[2]
    along<-fitrot$along[2]
    rho<-fitrot$rho[2]
    H11.2<-solve(invxtx[4:6,4:6])
    ahat<-fitrot$ahat13}
  if ((indside[1]==3)&(indside[2]==1)) {alat<- -fitrot$alat[2]
    along<-fitrot$along[2]+180
    rho<-fitrot$rho[2]
    H11.2<-fitrot$ahat13%*%solve(invxtx[4:6,4:6])%*%t(fitrot$ahat13)
    ahat<-t(fitrot$ahat13)}
  if ((indside[1]==2)&(indside[2]==3)) {alat<-fitrot$alat[3]
    along<-fitrot$along[3]
    rho<-fitrot$rho[3]
    H11.2<-invxtx[1:3,1:3]+invxtx[4:6,4:6]-invxtx[1:3,4:6]-invxtx[4:6,1:3]
    H11.2<-fitrot$ahat12%*%H11.2%*%t(fitrot$ahat12)
    H11.2<-solve(H11.2)
    ahat<-fitrot$ahat23}
  if ((indside[1]==3)&(indside[2]==2)) {alat<- -fitrot$alat[3]
    along<-fitrot$along[3]+180
    rho<-fitrot$rho[3]
    H11.2<-invxtx[1:3,1:3]+invxtx[4:6,4:6]-invxtx[1:3,4:6]-invxtx[4:6,1:3]
    H11.2<-fitrot$ahat13%*%H11.2%*%t(fitrot$ahat13)
    H11.2<-solve(H11.2)
    ahat<-t(fitrot$ahat23)}
} else {
#simple spreading ridge
  if ((indside[1]==1)&(indside[2]==2)) {alat<-fitrot$alat
    along<-fitrot$along
    rho<-fitrot$rho
    H11.2<-fitrot$H11.2
    ahat<-fitrot$ahat}
  if ((indside[1]==2)&(indside[2]==1)) {alat<- -fitrot$alat[1]
    along<-fitrot$along[1]+180
    rho<-fitrot$rho[1]
    H11.2<-fitrot$ahat%*%fitrot$H11.2%*%t(fitrot$ahat)
    ahat<-t(fitrot$ahat)}
}

rotmat<-ROTGEN2(plat,plong,glat,glong)
uxx<-as.double(1:3)
uxx<-.Fortran("trans1",alat,along,uxx)[[3]]
uxx<-rotmat%*%uxx
temp<-.Fortran("trans6",uxx,alat,along)
alat<-temp[[2]]
along<-temp[[3]]
ahat<-rotmat%*%ahat
H11.2<-rotmat%*%H11.2%*%t(rotmat)

qhat<-.Fortran("trans3",alat,along,rho,qhat)[[4]]
ahatm<-matrix(,nrow=3,ncol=4)
      ahatm[1,1]<--qhat[2]
      ahatm[1,2]<-qhat[1]
      ahatm[1,3]<-qhat[4]
      ahatm[1,4]<--qhat[3]
      ahatm[2,1]<--qhat[3]
      ahatm[2,2]<--qhat[4]
      ahatm[2,3]<-qhat[1]
      ahatm[2,4]<-qhat[2]
      ahatm[3,1]<--qhat[4]
      ahatm[3,2]<-qhat[3]
      ahatm[3,3]<--qhat[2]
      ahatm[3,4]<-qhat[1]
qbing<-4*t(ahatm)%*%H11.2%*%ahatm
qbing<-CONVMAT2(4,qbing)
lhat<-as.double(0)
if (is.finite(fitrot$df)) fk<-qf(alpha,df1,fitrot$df)*df1 else fk<-qchisq(alpha,df1)
fk<-fk/fitrot$kappahat
kmax<-as.integer(400)
workbd<-matrix(as.double(0),nrow=kmax,ncol=6)
nlong<-as.integer(201)
nlat<-as.integer(101)
icase<-as.integer(0)
ncurve<-as.integer(0)
nptbd<-as.integer(c(0,0,0))
uppmat<-matrix(as.double(0),nrow=nlong,ncol=nlat)
lowmat<-uppmat
jer<-as.integer(0)
temp<-.Fortran("bingham",qbing,lhat,fk,uxx,kmax,workbd,icase,ncurve,nptbd,nlat,nlong,uppmat,lowmat,jer)
ncurve<-temp[[8]]
if (ncurve==0) {workbd<-NULL
nptbd<-0} else {nptbd<-temp[[9]][1:ncurve]
workbd<-temp[[6]][1:max(nptbd),1:(2*ncurve)]}
icase<-temp[[7]]
uppmat<-temp[[12]]
uppmat<-ifelse(uppmat<(-1),NA,uppmat)
if (icase==7) {lowmat<-NULL} else {
lowmat<-temp[[13]]
lowmat<-ifelse(lowmat<(-1),NA,lowmat)}

list(filename=fitrot$filename,names=pnames,indside=indside,alpha=alpha,
kappahat=fitrot$kappahat,df=fitrot$df,alat=alat,along=along,rho=rho,
ahat=ahat,H11.2=H11.2,icase=icase,ncurve=ncurve,nptbd=nptbd,
bound=workbd,uppmat=uppmat,lowmat=lowmat,rotmat=rotmat,ier=temp[[14]],call.b=call,date.b=date())}