//calc2.js Bob Hanson adapted from Hans Reich  Laoccon Program

//CalcSpectrum(Nuclei,MinInt,mShifts,mJcoupl,gFreq,gInten)

Nuclei=0      //number of nuclei
MinInt=0.01   //the smallest peak height to handle
vs1=4.5
isdetail=false

gMAXLINE=5000      //maximum number of lines in a spectrum
gMAXNUCLEI=7

//Below the arrays for Full Spectrum Calculations
//note: array[0] will not be used

mJcoupl=new Array(gMAXNUCLEI+1)  //couplings in hz
for(var i=1;i<gMAXNUCLEI;i++)mJcoupl[i]=new Array()
mShifts=new Array(gMAXNUCLEI)    //shifts in hz


//outputs
gInten= new Array()
gFreq=new Array()
gWhich=new Array()  //assigment
NFreq=0 //number of lines

function CalcSpectrum(Nuclei,MinInt,mShifts,mJcoupl,gFreq,gInten){  //Seven Spin Calculation
	//---------------First fill the mShifts array (chemical shift in Hz)
	//[m][n], m<n (top right) half of coupling array is for +-J
	//[m][n], m>n (bottom left) half is for  +-J/2 just for speed

	for(var i=1;i<=Nuclei;i++){
		mJcoupl[i][i]=0
		for(var j=i+1;j<=Nuclei;j++)mJcoupl[j][i]=mJcoupl[i][j]*0.5
	}
	//-----------------Next calculate the stick spectrum
	//Lines outside the freq range are discarded
	//remember, freq increases right to left
	//This will fill the arrays gInten() and gFreq() to produce a stick spectrum
	//there are NFreq lines
	return CalcLineSecOrder()
}


function Sqr(x){return Math.pow(x,0.5)}

function Abs(x){return Math.abs(x)}

twotothe=new Array(1,2)
for (var i=2;i<16;i++)twotothe[i]=twotothe[i-1]*2

function ispowerof2(x){
	var p=0
	for(var m=0;m<16;m++){
		p=twotothe[m]
		if(p>x)return -1
		if(x==p)return m
	}
	return -1
}

function isUP(x,m){return (x&twotothe[m]?1:0)}

function nUP(x){
	var n=0
	for(var m=0;m<16;m++)n+=isUP(x,m)
	return n
}

function CalcLineSecOrder(){
	var i=0
	var ipt=0
	var n=0
	var J=0
	var JP=0
	var jd=0
	var m=0
	var k=0
	var n=0
	var DatPoint=0
	var na=0
	var nb=0
	var moa=0
	var mob=0
	var JL=0
	var kl=0
	var Ma=0
	var Mb=0
	var Ja=0
	var JB=0
	var npb=0
	var npa=0
	var JC=0
	var ka=0
	var kb=0
	var kd=0
	var FR=0
	var isum=0
	var sout=""
	var ptset=new Array()
	var No=new Array()
	var spins=new Array()
	var spinfac=new Array()
	var x1=new Array(36)
	for(var i=1;i<36;i++)x1[i]=new Array()
	var Vb=new Array(36)
	for(var i=1;i<36;i++)Vb[i]=new Array()
	var Va=new Array(36)
	for(var i=1;i<36;i++)Va[i]=new Array()
	
	var TempDat=new Array()       //Holds all Vb arrays
	var e=new Array()             //energies
	
	var nsqr=Math.pow(2,Nuclei)
	
	if(isdetail)sout+="Shifts, couplings:\n"+showmatrix(mShifts,0)+showmatrix(mJcoupl,Nuclei,1)
	
	//set up binomial series count No[]=(1,3,3,1) eg.
	//ptset points the the first of a set
	
	No[1]=1
	ptset[1]=1
	for(J=1;J<=Nuclei;J++){
		JP=J+1
		jd=Nuclei+1-J
		No[JP]=No[J]*jd/J
		ptset[JP]=ptset[J]+No[J]
	}
	
	//set up spins in sets of how many nuclei ar UP
	//the first set has none up, second 2 up, third 3 up, etc.
	//couplings will only be considered then between adjacent sets
	//as those have only one spin change
	
	spins[1]=0
	for(J=1;J<nsqr;J++){
		i=nUP(J)+1
		spins[ptset[i]]=J
		ptset[i]++
	}
	for(J=2;J<=Nuclei+1;J++)ptset[J]-=No[J]
	
	if(isdetail)sout+="SPINS:\n"+showmatrix(spins,0)
	if(isdetail)sout+="SPIN Sets:\n"+showmatrix(ptset,0)
	
	//set up factor matrix indicating +-1/2 based on each spin possibility
	
	for(k=1;k<=nsqr;k++){
		n=spins[k]
		spinfac[k]=new Array()
		for(m=0;m<Nuclei;m++)spinfac[k][m+1]=0.5-isUP(n,m)
	}
	
	if(isdetail)sout+="SPIN FACTOR MATRIX:\n"+showmatrix(spinfac,nsqr,1)
	
	//set up submatrices and diagonalize them
	
	DatPoint=0
	na=1
	while (na-Nuclei<=1){
		moa=No[na]
		
		if(isdetail)sout+="na, No[na]:"+na+" " + moa+"\n"
		
		for(J=1;J<=moa;J++){
			JL=ptset[na]+J-1
			for(k=J;k<=moa;k++){
				kl=ptset[na]+k-1
				if (J==k){
					x1[J][J]=0
					for(m=1;m<=Nuclei;m++)x1[J][J]+=spinfac[JL][m]*mShifts[m]
					for(m=1;m<Nuclei;m++){
						for(n=m+1;n<=Nuclei;n++)x1[J][J]+=spinfac[JL][m]*spinfac[JL][n]*mJcoupl[m][n]
					}
				}else{
					Ma=0;Mb=0
					for(m=1;m<=Nuclei;m++){
						if (spinfac[JL][m]*spinfac[kl][m]<0){
							if(Mb){Mb=0;break}
							if(Ma){Mb=m}else{Ma=m}
						}
					}
					
					if(isdetail)sout+="JL kl "+JL+" "+kl+" "+Ma+" "+Mb+" "+spinfac[JL][m]+" "+spinfac[kl][m]+"\n"
					
					x1[J][k]=x1[k][J]=(Mb?mJcoupl[Mb][Ma]:0)
				}
			}
		}
		
		if(isdetail)sout+="SUBMATRIX X1 in:\n"+showmatrix(x1,moa,1)
		
		if(!DiagMatrix(x1,Vb,moa))return false
		
		if(isdetail)sout+="SUBMATRIX X1 out:\n"+showmatrix(x1,moa,1)
		if(isdetail)sout+="SUBMATRIX Vb out:\n"+showmatrix(Vb,moa,1)
		
		for(J=1;J<=moa;J++){
			JL=J+ptset[na]-1
			e[JL]=x1[J][J]
		}
		
		if(isdetail)sout+="Energies e:\n"+showmatrix(e,0)
		
		for(J=1;J<=moa;J++){for(k=1;k<=moa;k++){
				TempDat[DatPoint]=Vb[J][k]
				DatPoint++
		}}
		na++
	}
	
	// get all frequencies
	
	NFreq=0
	DatPoint=0
	Vb[1][1]=TempDat[DatPoint]
	DatPoint++
	na=0
	while(Nuclei>na){
		na++
		nb=na+1
		moa=No[na]
		mob=No[nb]
		for(JB=1;JB<=mob;JB++){
			npb=JB+ptset[nb]-1
			for(Ja=1;Ja<=moa;Ja++)Va[Ja][JB]=0
			for(Ja=1;Ja<=moa;Ja++){
				npa=Ja+ptset[na]-1
				//must differ by one spin only
				if(ispowerof2(spins[npb]-spins[npa])>=0){
					
					if(isdetail)sout+="\n linked: "+spins[npb]+" "+spins[npa]
					
					for(JC=1;JC<=moa;JC++)Va[JC][JB]+=Vb[Ja][JC]
				}
			}
		}
		//load next Vb matrix
		for(J=1;J<=mob;J++){for(k=1;k<=mob;k++){
				Vb[J][k]=TempDat[DatPoint]
				DatPoint++
		}}
		
		if(isdetail)sout+="Freq Matrices: Va,Vb:\n"+showmatrix(Va,moa,1)+showmatrix(Vb,mob,1)
		
		for(ka=1;ka<=moa;ka++){
			npa=ka+ptset[na]-1
			for(kb=1;kb<=mob;kb++){
				npb=kb+ptset[nb]-1
				i=0
				for(kd=1;kd<=mob;kd++)i+=Va[ka][kd]*Vb[kd][kb]
				i*=i
				if (i>=MinInt){
					NFreq++
					if (NFreq>gMAXLINE)return false
					isum+=i
					gFreq[NFreq]=e[npa]-e[npb]
					
					if(isdetail)sout+="Frequency npa npb e" + NFreq+": "+gFreq[NFreq]+" "+npa+" "+npb+" "+e[npa]+" "+e[npb] + " "+i+"\n"
					
					gInten[NFreq]=i
					gWhich[NFreq]=ispowerof2(spins[npb]-spins[npa])
				}
			}
		}
		
		if(isdetail)sout+="gFreq,gInten,gWhich:"+NFreq+"\n"+showmatrix(gFreq,0)+showmatrix(gInten,0)+showmatrix(gWhich,0)
		
	}
	
	if(isdetail)document.info.list.value=sout
	
	//normalize
	k=(isum?Nuclei*vs1/isum:0)
	for(i=1;i<=NFreq;i++)gInten[i]=gInten[i]*k
	return true
}

function DiagMatrix(M0, M1, dim){
	var nsqp=0
	var iterm=0
	var ipi=0
	var nmi=0
	var mati=0
	var MATJ=0
	var matk=0
	var matl=0
	var JM1=0
	var BIGX=0
	var BIGR=0
	var DEL=0
	var ROT=0
	var TS=0
	var r=0
	var AQ=0
	var BQ=0
	var F1=0
	var d=0
	nsqp=Math.floor(3*dim*dim/2)
	iterm=0
	nmi=dim-1
	
	//if(isdetail)alert(showmatrix(M0,dim,1))
	
	if (dim==0)return false
	if(dim==1){
		M1[1][1]=1
		return true
	}
	for(mati=1;mati<=dim;mati++){
		for(MATJ=1;MATJ<=dim;MATJ++)M1[mati][MATJ]=0
	}
	for(matk=1;matk<=dim;matk++) M1[matk][matk]=1
	//Find largest off-diagonal element
	while(nsqp>=iterm){
		iterm++
		BIGX=0
		for(MATJ=2;MATJ<=dim;MATJ++){
			JM1=MATJ-1
			for(mati=1;mati<=JM1;mati++){
				if (BIGX<Abs(M0[mati][MATJ])){
					BIGX=Abs(M0[mati][MATJ])
					matk=mati
					matl=MATJ
				}
			}
		}
		if (BIGX <= 0.000009999999){
			BIGR=0.00004
			for(mati=1;mati<=nmi;mati++){
				ipi=mati+1
				for(MATJ=ipi;MATJ<=dim;MATJ++){
					if (M0[mati][MATJ]){
						DEL=M0[mati][mati]-M0[MATJ][MATJ]
						if(DEL){
							ROT=Abs(M0[mati][MATJ]/DEL)
							if(ROT>BIGR){
								BIGR=ROT
								matk=mati
								matl=MATJ
							}
						}
					}
				}
			}
			//degenerate energy levels
			
			//if(isdetail)alert(BIGR+showmatrix(M0,dim,1))
			
			if (BIGR<=0.0000005)return true
		}
		//Compute Rotation
		TS=M0[matk][matl]*M0[matk][matl]
		DEL=M0[matk][matk]-M0[matl][matl]
		r=Sqr(Abs((DEL*DEL)+(4*TS)))
		AQ=Sqr(Abs((r+DEL)/(2*r)))
		if(AQ<0.707){
			BQ=-AQ
			AQ=Sqr(1-BQ*BQ)
		}else{
			BQ=-Sqr(1-AQ*AQ)
		}
		if (DEL/M0[matk][matl]<0)BQ=-BQ
		for(MATJ=1;MATJ<=dim;MATJ++){
			F1=M1[MATJ][matk]
			M1[MATJ][matk]=AQ*F1-BQ*M1[MATJ][matl]
			M1[MATJ][matl]=BQ*F1+AQ*M1[MATJ][matl]
			//Orthogonal rotation of matrix
			if(MATJ != matk && MATJ != matl){  //if either is equal, then don//t
				M0[matk][MATJ]=AQ*M0[MATJ][matk]-BQ*M0[MATJ][matl]
				M0[matl][MATJ]=BQ*M0[MATJ][matk]+AQ*M0[MATJ][matl]
				M0[MATJ][matk]=M0[matk][MATJ]
				M0[MATJ][matl]=M0[matl][MATJ]
			}
		}
		d=M0[matk][matk]+M0[matl][matl]
		M0[matk][matk]=(AQ*AQ)*M0[matk][matk]+(BQ*BQ)*M0[matl][matl]-2*AQ*BQ*M0[matk][matl]
		M0[matl][matl]=d-M0[matk][matk]
		M0[matl][matk]=0
		M0[matk][matl]=0
	}
	return true
}