#pragma rtGlobals=1		// Use modern global access method.

// Put the line
//            #include "Scalc"
//  at the top of your procedure window (but remove the leading "//" and spaces)
//  This will automaitcally include it in your experiment and make it available to you.
//  Or, you can just go to the "File" menu, and go down to the the "Open File" sub menu, and select "Procedure..."


Function Scalc(qkF,ene,Gc,XoRe,XoIm,XoEne)		// calculates S(q,w,G), given three waves representing Xo (real,imag,energy).
	Variable qkF		// q (kF) for Aluminum (kF = 1.75)
	Variable ene		// energy loss where to calculate
	Variable/C Gc		// complex G to use
	Wave XoRe			// wave with Re(Xo)
	Wave XoIm			// wave with Im(Xo)
	Wave XoEne			// wave with ĈE that goes with Xo, energy loss is positive

	Variable/C Xo		// complex value of Xo at ene
	Variable/C one = cmplx(1.,0.)
	Variable kF=1.75	// kF for Aluminum
	Variable vq			// v in (1/eV) assuming q in (kF)

	vq = 102.681 / (qkF*kF)^2	// compute the v (Coulomb potential)
	Xo = cmplx(interp(ene, XoEne, XoRe),interp(ene, XoEne, XoIm))
	return -imag(Xo / ( one - vq*( one-Gc)*Xo))
End


Macro Scalc_Example()		// An example it uses a polynomial approx for Xo at 1.7 kF for testing,
	Silent 1					// 		not needed for operation of Scalc, just an example

	//   Ignore the following 7 lines
	//
	Make/N=101/O Xr_test,Xi_test,Xe_test		// create a dummy Xo just for testing, normally you get this from Fleszar
	Xe_test = p/(numpnts(Xe_test)-1) * 100	// the energy points
	Make/O W_coef__={-0.028213,-0.00013953,0.0001222,-6.133e-06,1.5202e-07,-2.0403e-09,1.396e-11,-3.8024e-14}
	Xr_test= poly(W_coef__,Xe_test[p])
	Make/O W_coef__={0.00054337,-0.0021016,9.2358e-05,-2.0945e-06,2.7596e-08,-1.666e-10,5.1321e-14,2.3571e-15}
	Xi_test= poly(W_coef__,Xe_test[p])
	KillWaves/Z W_coef__


	// a simple call to Scalc
	//
	print "S(1.5kF,25eV,G=(1.5,0.1)) =",Scalc(1.5,25,cmplx(1.5,0.1),Xr_test,Xi_test,Xe_test)	// print S(q,w) at one point


	// creation of S(q,w,G=const) using Scalc
	//
	Make/N=21/O Stest1
	SetScale/I x 0,80,"eV", Stest1
	Stest1 = Scalc(2.1,x,cmplx(0.5,0.1),Xr_test,Xi_test,Xe_test)	// fill Stest1 for q=2.1kF


	// creation of S(q,w,G(E)) using Scalc
	//
	Make/N=21/O Stest2
	Make/N=21/O/C Gtest2
	SetScale/I x 0,80,"eV", Stest2, Gtest2
	Gtest2 = cmplx(1-0.012*x+0.00012*x^2,x/200)				// fill G wave with form of G
	Stest2 = Scalc(2.1,x,Gtest2(x),Xr_test,Xi_test,Xe_test)		// fill Stest2 for q=2.1kF

	display Stest1,Stest2											// and display the results
	append/R Gtest2
	ModifyGraph lsize=2,mode(Gtest2)=3,marker=19,msize=2,rgb(Stest1)=(0,0,65535)
	SetAxis/A/E=1
	Legend/N=text0/S=3/M
EndMacro