#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3

// calculate phase diagram using Kd paramaterization a la Langmuir and 
// Hanson (1981) in Advances in Physical Geochemistry Vol 1, 
// Thermodynamics of Minerals and Melt

// ---------- make a ternary grid

function makeGrid(divisions, strA, strB, strC)
	variable divisions // 100 divisions for 1% increments
	string strA, strB, strC
	
	variable Pnts=(divisions+1)*(divisions/2+1)
	make /o/n=(Pnts) $strA /Wave=componentA
	make /o/n=(Pnts) $strB /Wave=componentB
	make /o/n=(Pnts) $strC /Wave=componentC
	
	variable A,B,i=0
	for(A=0;A<divisions+1;A+=1)
		for(B=0;B<(divisions-A+1);B+=1)
			componentA[i]=A/divisions
			componentB[i]=B/divisions
			componentC[i]=1-componentA[i]-componentB[i]
			i+=1
		endfor
	endfor
	
// rounding is useful if you want nice decimal divisions
//	componentA=abs(RoundToDecimal(ceil(1+log(divisions)), componentA))
//	componentB=abs(RoundToDecimal(ceil(1+log(divisions)), componentB))
//	componentC=abs(RoundToDecimal(ceil(1+log(divisions)), componentC))
end
			
function RoundToDecimal(i, val)
	variable i, val
	
	return round(val*10^i)/10^i
end

function TernaryLiquidus(XA,XB)
	variable XA, XB
	
	Variable XC=1-XA-XB
	variable TLiqA, TLiqB, TLiqC
		
	// w[0] is the component
	// make /free w={1, XFo, XDi, XAn, 1890, 1392, 1553, 1387, 1274, 1327, 0.9, 0.42, 0.29}
	
	// wave w contains the component number, component proportions, 
	// component melting temperatures, eutectic temperatures, and eutectic
	// component proportions
	// w={component, XA, XB, XC, TmA, TmB, TmC, TeAB, TeBC, TeCA, XBeAB, XCeBC, XAeCA}
	
	make /free w={1, XA, XB, XC, 500, 450, 470, 250, 300, 275, 0.6, 0.4, 0.35}
	w[4,9]+=273.15 // C -> K
	
	variable LowLimit=100, HighLimit=5000
	
	if (XA<1e-5)
		TLiqA=-inf
	elseif(XA>1-1e-5)
		return w[4]-273.15
	else
		w[0]=0
		findroots /Q/H=(HighLimit)/L=(LowLimit)/Z=(XA) TernaryX, w
		TLiqA=v_root-273.15
	endif
		
	if (XB<1e-5)
		TLiqB=-inf
	elseif(XB>1-1e-5)
		return w[5]-273.15
	else
		w[0]=1
		findroots /Q/H=(HighLimit)/L=(LowLimit)/Z=(XB) TernaryX, w
		TLiqB=v_root-273.15
	endif
	
	if (XC<1e-5)
		TLiqC=-inf
	elseif(XC>1-1e-5)
		return w[6]-273.15
	else
		w[0]=2
		findroots /Q/H=(HighLimit)/L=(LowLimit)/Z=(XC) TernaryX, w
		TLiqC=v_root-273.15
	endif
		
	return max(TLiqA,TLiqB,TLiqC)
end

function TernaryX(w, T)
	wave w
	variable T
			
	// caculate proportion of component w[0] at temperature T
	// in ternary eutectic defined by parameters in w
	// assuming ideal mixing of ternary components in liquid
	// no solid solution between component phases
		
	// assume ternary components 0, 1 and 2 are arranged in clockwise fashion
	
	// use an Arrhenius relationship
	// log Kd = alpha/T + beta
	
	// alpha1 and beta1 for mixing with the component on the clockwise side,
	// alpha2 and beta2 counterclockwise.
	
	// mix linearly in logKd - 1/T space	
	
	variable alpha1, beta1, alpha2, beta2, alphaMix, betaMix
	variable Tm, T_eutectic1, T_eutectic2, X_eutectic1, X_eutectic2
	variable X1, X2 // proportions of components at the other end of join 1 and join 2
	variable component=w[0]
	
	Tm=w[4+component]
	
	// clockwise	side of ternary
	T_eutectic1=w[7+component]
	// X(component) in eutectic composition
	X_eutectic1=1-w[10+component]	
	X1=w[1+component+1-3*(component==2)]
	alpha1=-log(X_eutectic1)/(1/T_eutectic1-1/Tm)
	beta1=-alpha1*1/Tm
	
	// anticlockwise side of ternary
	T_eutectic2=w[7+component-1+3*(component==0)]
	X_eutectic2=w[10+component-1+3*(component==0)]
	X2=w[1+component-1+3*(component==0)]
	alpha2=-log(X_eutectic2)/(1/T_eutectic2-1/Tm)
	beta2=-alpha2*1/Tm
	
	alphaMix=X1/(X1+X2)*alpha1+X2/(X1+X2)*alpha2
	betaMix=X1/(X1+X2)*beta1+X2/(X1+X2)*beta2
	
	return 1/10^(betaMix+alphaMix/T)
end


// here's how we do solid solutions (divariant melting loops)

function AlbiteAnorthiteLiquidus(An)
	variable An
	
	make /free w={0}
	
	findroots /Q/H=1830/L=1350/Z=(An) AnAbLiquidXAn, w
	return v_root-273.15
end

function AlbiteAnorthiteSolidus(An)
	variable An
	
	make /free w={0}
	
	findroots /Q/H=1830/L=1350/Z=(An) AnAbSolidXAn, w
	return v_root-273.15
end

function AnAbLiquidXAn(w, T)
	wave w; variable T
	
	variable KdAb=10^(2512/T-1.85)
	variable KdAn=10^(6358/T-3.49)
	return (1-KdAb)/(KdAn-KdAb)
end

function AnAbSolidXAn(w, T)
	wave w 
	variable T
	
	// first find liq composition at T
	variable liqAn=AnAbLiquidXAn(w, T)
	// now find solid composition
	variable KdAn=10^(6358/T-3.49)
	return liqAn*KdAn	
end