vibscale=.8  		//default magnitude
vectoronlyscale=20	//factor when vectors only

//vmag=1      		//user-selectable magnitude
//isvectorsonly=0 	//set to 1 if vectors only



defaultvibscript="wireframe .13;spacefill 25%;background white;frame 1;vectors on;color vectors cyan;vibration on;set frank off"

function createjmoldata(Structure,q){
 var format="| nnn        nnn         nnnnnnnnnn  nnnnnnnnnn  nnnnnnnnnn\n"
 var sline=""
 var s="| Entering Gaussian System\n"
 s+="| This is part of the Gaussian 94(TM) system of programs. \n"
 s+="| Standard orientation:\n"
 s+="|\n"
 s+="|\n"
 s+="|\n"
 s+="|\n"
 
 for(var i=0;i<Structure.Atno.length;i++){
	sline=format.replace(/nnn/,flushRight(i+1,3)).replace(/nnn/,flushRight(Structure.Atno[i],3))
	sline=sline.replace(/nnnnnnnnnn/,flushRight(roundoff(Structure.Coords[i][0],6),10))
	.replace(/nnnnnnnnnn/,flushRight(roundoff(Structure.Coords[i][1],6),10))
	.replace(/nnnnnnnnnn/,flushRight(roundoff(Structure.Coords[i][2],6),10))
	s+=sline
 }
 
 s+="| ------------\n"

 format="| nnn       nnnnnn nnnnnn nnnnnn\n"
        //   1        -0.12   0.00   0.00  

 s+="| Harmonic frequencies\n"
 s+="| Frequencies -- \n"
 s+="|\n"
 s+="|\n"
 s+="|\n"
 s+="|\n"
 s+="|\n"
 s+="|\n"

 for(var i=0;i<Structure.Atno.length;i++){
	sline=format.replace(/nnn/,flushRight(i+1,3))
	sline=sline.replace(/nnnnnn/,flushRight(roundoff(Structure.Vectors[q][i][0],2),6))
	.replace(/nnnnnn/,flushRight(roundoff(Structure.Vectors[q][i][1],2),6))
	.replace(/nnnnnn/,flushRight(roundoff(Structure.Vectors[q][i][2],2),6))
	s+=sline
 }
 return s
}

function createjmoldata(Structure,q){ 
var s="|" + Structure.Atno.length + "\n|CoolMolecules Vibration Generator\n"
 for(var i=0;i<Structure.Atno.length;i++){
	s += "|" + Elements[Structure.Atno[i]].sym + " " 
	+ Structure.Coords[i][0] + " " + Structure.Coords[i][1] + " " + Structure.Coords[i][2]+ " " 
	+ Structure.Vectors[q][i][0] + " " + Structure.Vectors[q][i][1] + " " + Structure.Vectors[q][i][2] + "\n" 
}
 
return s
}


function getvibration(modelname,modelinfo,modename,isvectorsonly){
 var scalefactor=(isvectorsonly?20:1)*vmag*vibscale
 Model=Models[modelname]  //global
 var Structure=new Array()
 Structure.Coords=new Array()
 Structure.Vectors=new Array()
 Structure.Atno=new Array()
 readinfo(modelinfo)
 setcoord(Structure)

 dim=Model.Modes[modename].Vectors.length
 var l=Structure.Coords.length

 for(var i=0;i<l;i++)for(var j=0;j<3;j++)if(Math.abs(Structure.Coords[i][j])<1E-5)Structure.Coords[i][j]=0

 for(var i=0;i<dim;i++)setvectors(Structure,Model.Modes[modename],i,scalefactor)

 for(var i=0;i<Structure.Vectors[0].length;i++)Structure.Vectors[0][i].moves=false
 for(var i=1;i<Structure.Vectors[0].length;i++)for(var j=0;j<dim;j++)if(Model.Modes[modename].Vectors[j][i][1]!=0)Structure.Vectors[0][i].moves=true

 for(var i=0;i<dim;i++)angular(Structure,Model.Modes[modename],i) 

 for(var i=0;i<dim;i++)vcentermass(Structure,Model.Modes[modename],i,dim)
 for(var i=0;i<dim;i++)for(var j=0;j<l;j++)for(var k=0;k<3;k++)if(Math.abs(Structure.Vectors[i][j][k])<1E-5)Structure.Vectors[i][j][k]=0

 var lmax=0
 var maxlength=(isvectorsonly?1.6:.08)*vmag
 for(var i=0;i<dim;i++)for(var j=0;j<Structure.Vectors[0].length;j++)if(vlength(Structure.Vectors[i][j])>lmax)lmax=vlength(Structure.Vectors[i][j])
 for(var i=0;i<dim;i++)for(var j=0;j<Structure.Vectors[0].length;j++)copyarray(vscale(Structure.Vectors[i][j],maxlength/lmax),Structure.Vectors[i][j])
 var D=new Array()
 D[0]=new Array()
 for(var i=0;i<dim;i++)for(var j=0;j<l;j++){
	D[0][D[0].length]=[0,0,0]
	copyarray(Structure.Vectors[i][j],D[0][D[0].length-1])
 }
 for(var i=0;i<D[0].length;i++){
	Structure.Vectors[0][i]=[0,0,0]
	copyarray(D[0][i],Structure.Vectors[0][i])
 }
 displacement=Model["X-A"]*.04
 for(var i=1;i<dim;i++)for(var j=0;j<l;j++)Structure.Atno[Structure.Atno.length]=Structure.Atno[j]
 var C=new Array()
 for(var i=0;i<l;i++)C[i]=[0,0,0]
 for(var i=0;i<l;i++)copyarray(Structure.Coords[i],C[i])
 for(var i=1;i<dim;i++){
	for(var j=0;j<l;j++){
			C[C.length]=[0,0,0]
			copyarray(Structure.Coords[j],C[C.length-1])
			C[C.length-1][0]+=(displacement*i)
			
	}
 }
 for(var i=0;i<C.length;i++){
	Structure.Coords[i]=[0,0,0]
	copyarray(C[i],Structure.Coords[i])
 }
 Model.Structure=Structure

 var dim1=1
 Jmolcode=new Array()
 for(var i=0;i<dim1;i++)Jmolcode[i]=createjmoldata(Structure,i)
 return Jmolcode[0]
}

function readinfo(modelinfo){
 
 var formula=""
 var atomcount=0
 var lastatom=""
 var atom=""
 var Modeldata=modelinfo.split("?")
 for(var i=0;i<Modeldata.length;i++)Modeldata[i]=Modeldata[i].split("=")
 var S=new Array()

 for(var i=0;i<Model.XAB.length;i++){
	atom=Model.XAB[i]
	for(var j=0;j<Modeldata.length;j++){
		if(Modeldata[j][0]==atom)Model.Atoms[i]=Modeldata[j][1]
	}
 }

 for(var i=0;i<Model.Atoms.length;i++){
	if(Model.Atoms[i]!==lastatom){
		if(i!=0&&atomcount>1)formula+=atomcount
		formula+=Model.Atoms[i]
		atomcount=1
	}else{
		atomcount++
	}
	lastatom=Model.Atoms[i]
 }
 if(atomcount>1)formula+=atomcount
 Model.formula=formula
 for(var i=0;i<Modeldata.length;i++){
	if(Modeldata[i][0]=="d")Model[Modeldata[i][1]]=parseFloat(Modeldata[i][2])
 }
}

function setcoord(Structure){
 var rotxyz="x"
 var angle=0
 var dist=0
 var A=[0,0,0]
 var S=new Array()
 Structure.Coords[0]=[0,0,0]
 Structure.Atno[0]=AtomNum[Model.Atoms[0]]
 for(var i=1;i<Model.Coords.length;i++){
	Structure.Atno[i]=AtomNum[Model.Atoms[i]]
	Structure.Coords[i]=new Array()
	dist=Model[Model.Coords[i][0]]
	A=[dist,0,0]
	if(A[0]>50)A[0]/=100
	for(var j=1;j<Model.Coords[i].length;j++){
		S=Model.Coords[i][j].split(":")
		rotxyz=S[0]
		angle=S[1].replace(/\]/g,"']").replace(/\[/g,"Model['")
		angle=eval(angle)
		if(rotxyz=="z"){
			A=vrotatez(A,angle)
		}else if(rotxyz=="y"){
			A=vrotatey(A,angle)
		}else{
			A=vrotatex(A,angle)
		}		
	}
	copyarray(A,Structure.Coords[i])
 }
}

function setvectors(Structure,Mode,q,scalefactor){
 var rotxyz="x"
 var angle=0
 var dist=0
 var factor=.1
 var i0=0
 var A=[0,0,0]
 var S=new Array()
 Structure.Vectors[q]=new Array()
 Structure.Vectors[q][0]=[0,0,0]
 for(var i=1;i<Mode.Vectors[q].length;i++){
	Structure.Vectors[q][i]=new Array()
	i0=Mode.Vectors[q][i][0]
	if(i0=="l")copyarray([0,0,0],Structure.Vectors[q][i])
	if(i0=="a"){
		var Proj=new Array()
		var dot=0
		var H=new Array()
		var K=new Array()	
		var L=[0,0,0]
		for(var j=1;j<Mode.Vectors[q][i].length-1;j+=2){
			copyarray(Structure.Coords[Mode.Vectors[q][i][j]],K)
			K=vscale(K,Mode.Vectors[q][i][j+1])
			L=addarray(L,K)
		}
		copyarray(Structure.Coords[i],H)
		for(var j=0;j<3;j++)dot+=L[j]*H[j]
		Proj=vscale(H,dot/(vlength(H)*vlength(H)))
		K=addarray(L,vscale(Proj,-1))
		K=vscale(K,scalefactor*Mode.Vectors[q][i][Mode.Vectors[q][i].length-1]/vlength(K))
		copyarray(K,Structure.Vectors[q][i])
	}
	if(typeof(i0)=="number"){

		factor=Mode.Vectors[q][i][1]
		copyarray(Structure.Coords[i0],A)
		A=vscale(A,factor*scalefactor/vlength(A))
		for(var j=2;j<Mode.Vectors[q][i].length;j++){
			if(Mode.Vectors[q][i][j]!="l"){
				S=Mode.Vectors[q][i][j].split(":")
				rotxyz=S[0]
				angle=S[1].replace(/\]/g,"']").replace(/\[/g,"Model['")
				angle=eval(angle)
				if(rotxyz=="z"){
					A=vrotatez(A,angle)
				}else if(rotxyz=="y"){
					A=vrotatey(A,angle)
				}else{
					A=vrotatex(A,angle)
				}		
			}
		}
		copyarray(A,Structure.Vectors[q][i])
	}
 }
}

function vcentermass(Structure,Mode,q,dim){
 var Newcoords=new Array()
 for(var i=1;i<Structure.Coords.length;i++){
	Newcoords[i]=addarray(Structure.Coords[i],Structure.Vectors[q][i])
 }

 var Centersm=[0,0,0]
 var sumatmass=arraysum(Structure.Atmass)
 for(var j=0;j<3;j++){
	 for(var i=0;i<Structure.Atmass.length;i++){
 		Centersm[j]+=Structure.Atmass[i]*Structure.Coords[i][j] 
	 }
 	 Centersm[j]/=sumatmass
 }
 for(var j=0;j<3;j++){
	 var t=Structure.Atmass[0]
	 var s=Centersm[j]*sumatmass
	 for(var i=1;i<Structure.Vectors[q].length;i++){
		if(Structure.Vectors[0][i].moves){
			s-=Newcoords[i][j]*Structure.Atmass[i]
		}else{
			t+=Structure.Atmass[i]
			s-=Structure.Coords[i][j]*Structure.Atmass[i]
		}	
	 }
	 Structure.Vectors[q][0][j]=s/t
	 for(var i=1;i<Structure.Coords.length;i++){ 
 		if(!Structure.Vectors[0][i].moves)Structure.Vectors[q][i][j]=s/t
	 }
 }
}


function angular(Structure,Mode,q){
 Structure.Atmass=new Array()
 for(var i=0;i<Structure.Atno.length;i++)if(Structure.Atno[i])Structure.Atmass[i]=Elements[Structure.Atno[i]].atmass
 var Q=[0,0,0]
 var P=new Array()
 var R=new Array()
 var T=new Array()
 var L=[0,0,0]
 var numatoms=Structure.Atmass.length
 for(var i=1;i<numatoms;i++){
 	if(Mode.Vectors[q][i][Mode.Vectors[q][i].length-1]!="l")L=addarray(L,vscale(vcross(Structure.Coords[i],Structure.Vectors[q][i]),Structure.Atmass[i]))
 }
 for(var k=0;k<3;k++)if(Math.abs(L[k])<1E-5)L[k]=0
 var b=0
 for(var i=1;i<numatoms;i++)if(Mode.Vectors[q][i][0]=="l")b++
 for(var i=1;i<numatoms;i++){
	if(Mode.Vectors[q][i][Mode.Vectors[q][i].length-1]=="l"&&Mode.Vectors[q][i].length>1){
		P[i]=new Array()
		copyarray(Structure.Vectors[q][i],P[i])
	}
	if(Mode.Vectors[q][i][0]=="l"){
		R=vscale(Structure.Coords[i],Structure.Atmass[i])
		T[0]=[0,-R[2],R[1],-L[0]/b]
		T[1]=[-R[2],0,R[0],-L[1]/b]	
		T[2]=[-R[1],R[0],0,-L[2]/b]
		T=rref(T)
		for(var j=0;j<3;j++){
			if(T[j][j]==1){
				Structure.Vectors[q][i][j]=T[j][3]
			}else{
				Structure.Vectors[q][i][j]=0
			}
		}	
	}
 }
 var v=0
 for(var i=0;i<numatoms;i++)if(P[i])P[i]=vscale(P[i],Structure.Atmass[i])
 for(var i=0;i<numatoms;i++)if(P[i])Q=addarray(vcross(P[i],Structure.Coords[i]),Q)
 for(var i=0;i<3;i++)if(Math.abs(Q[i])>1E-7)v=L[i]/Q[i]
 for(var i=0;i<numatoms;i++)if(P[i])copyarray(vscale(Structure.Vectors[q][i],v),Structure.Vectors[q][i])
}





function rref(A){
 var m=A.length
 var D=new Array()
 for(var i=0;i<m;i++){
	if(A[i][i]==0){
		for(var j=i+1;j<m;j++){
			if(A[j][i]!=0){		
				copyarray(A[i],D)
				copyarray(A[j],A[i])
				copyarray(D,A[j])
				j=m
			}
		}
	}
	if(A[i][i]){
		var f=A[i][i]
		for(var j=0;j<A[i].length;j++)A[i][j]/=f
	}
	for(var k=0;k<m;k++)if(k!=i)copyarray(addarray(vscale(A[i],-A[k][i]),A[k]),A[k])
 }
 for(var i=0;i<m;i++)for(var j=0;j<A[0].length;j++)A[i][j]=Math.round(1E7*A[i][j])/1E7
 return(A)
}

function showmatrix(A){
 var s=""
 for(var i=0;i<A.length;i++)s+=A[i]+"\n"
 return(s)
}

function vcross(A,B){
 return([A[1]*B[2]-A[2]*B[1],A[0]*B[2]-A[2]*B[0],A[0]*B[1]-A[1]*B[0]])
}

function vlength(A){
 return Math.sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2])
}

function vrotatez(A,deg){
 var rad=deg*Math.PI/180
 return [A[0]*Math.cos(rad)-A[1]*Math.sin(rad),A[0]*Math.sin(rad)+A[1]*Math.cos(rad),A[2]]
}

function vrotatey(A,deg){
 var rad=deg*Math.PI/180
 return [A[0]*Math.cos(rad)+A[2]*Math.sin(rad),A[1],-A[0]*Math.sin(rad)+A[2]*Math.cos(rad)]
}

function vrotatex(A,deg){
 var rad=deg*Math.PI/180
 return [A[0],A[1]*Math.cos(rad)-A[2]*Math.sin(rad),A[1]*Math.sin(rad)+A[2]*Math.cos(rad)]
}

function vscale(A,f){
 return [A[0]*f,A[1]*f,A[2]*f]
}


//from code.js

function addarray(C,D){
 var E=new Array()
 for(var i=0;i<C.length;i++)E[i]=C[i]+D[i]
 return E
}
 
function arraysum(B){
 var s=0
 for(var i=0;i<B.length;i++)s+=B[i]
 return s
}

function copyarray(A,B){
 for(var i=0;i<A.length;i++)B[i]=A[i]
}

function flushRight(s,n){
	var sp="                                       "+s
	return sp.substring(sp.length-n,sp.length)
}

function roundoff(x,ndec){

	//ndec>0 for fixed; ndec<0 for scientific
	var s=""

	if(isNaN(x) || x==0)return "0"

	if(ndec==0)return ""+Math.round(x)

	var neg=(x<0?"-":"")

	var xs=Math.abs(x)+""

	var i=(xs.indexOf("E") & xs.indexOf("e"))

	if(ndec<0 && i<0){

		xs=roundoff(Math.abs(x)*1e-100,-ndec)

		i=(xs.indexOf("E") & xs.indexOf("e"))

		var e=parseInt(xs.substring(i+1,xs.length))+100

		s=neg+xs.substring(0,i)+(e!=0?"E"+e:"")

		return s

	}

	if (i>0){

		s=roundoff(xs.substring(0,i),Math.abs(ndec)-1)+"E"+xs.substring(i+1,xs.length)

		if(s.indexOf("10.")==0){

			i=(s.indexOf("E") & s.indexOf("e"))

			s="1"+s.substring(2,i+1)+(parseInt(s.substring(i+1,s.length))+1)

		}

		return neg+s

	}

	i=xs.indexOf(".")

	if(i<0){xs+=".";i+=xs.length}

	xs+="000000000"

	s="."+xs.substring(i+1+ndec,xs.length)

	xs=xs.substring(0,i)+xs.substring(i+1,i+1+ndec)

	var add1=(xs.charAt(0)=="0")

	if(add1)xs="1"+xs

	xs=parseInt(xs)+Math.round(parseFloat(s))+""

	if(add1)xs=xs.substring(1,xs.length)

	xs=xs.substring(0,xs.length-ndec)+"."+xs.substring(xs.length-ndec,xs.length)

	if(xs.indexOf(".")==0)xs="0"+xs

	return neg+xs

}