setCaseIgnored false

to nmr_carterspectra
	setUpdateGraph "false
	hideTurtle
	clearScreen
	WindowMode
	setpensize 0
;	disableTextureFont
	enablePolySmooth
	ch=0
nmr
stop
	forever 
	[
		clearScreen
		;pr repcount
		catch "error [ignore runResult [nmr]]
		err=Error
		if :err != []
		[	clearText
			print :err.2
			waitMS 200
		]
		updategraph
		updateVars
		dispatchMessages
		waitMS 100
		if ch==Char 27 [stop]
	]
end

to nmr
	ni=1000
	x=rSeqFA 400 -400 :ni

	spectrum=readfile "|../La2RuO5/56Mhz/0_700K_ausw.dat|
	specx=:spectrum.1
	specy=:spectrum.2
	specy=:specy/(max :specy)
	specn=count :specy
	specy_cor=:specy-rSeqFA 0.0 0.05 specn
	specy=:specy_cor*400-200

;	B=specx
;	nni=count specx
;	BI=Array nni
;	repeat nni
;	[	i=repcount
;		BI.i=list B.i specy_cor.i
;	]
;	foreach BI [pr ?]

	spectrum2=readfile "|../La2RuO5/64Mhz/1K_2_ausw.dat|
	specx2=:spectrum2.1/64*56
	specy2=:spectrum2.2
	specy2=:specy2/(max :specy2)*400-200

;	B=specx2
;	nni=count specx2
;	BI=Array nni
;	repeat nni
;	[	i=repcount
;		BI.i=list B.i spectrum2.(2).i
;	]
;	foreach BI [pr ?]
	
	for [p 0 0]
	[	intensity=FloatArray :ni

		I=7/2
		MHz=1e6
		Tesla=1
		gamma=6.061*:MHz/:Tesla
		markerx=56.25*:MHz
		nu0=56.25*:MHz
		nuref=56*:MHz
		nuQ=4.8*:MHz
		eta=0.85 
		expo=0
		K1=0 ;0.0045
		K2=0
		ntheta=200
		nphi=200
		I0=1
		Inu0=:nu0
		nnu=3*:Tesla*:gamma
		smoothwidth=0.225*:MHz
		compwidth=10*:MHz
		maxy=200
		maxx=400
		ticksx=5
		tick=15
		txtx=25

		xmax=(-:nu0+:nnu/2)/:gamma
		xmin=(-:nu0-:nnu/2)/:gamma
		w=:xmax-:xmin
		wq=(:xmax-:xmin)*2^(int -0.5+(ln :w/(:xmax-:xmin))/ln 2)
		markerx=:markerx/:gamma

		nmrloop
		if ch==char 27 [break]
	]
;	B=rSeqFA -xmin -xmax ni
;	BI=Array ni
;	repeat ni
;	[	i=repcount
;		BI.i=list B.i smoothedIntensity.i
;	]
;	foreach BI [pr ?]
end

to nmrloop
	repeat 1
	[
;		computeIntensityRatios
;stop
		intensity=NMRCarterIntensity 
			:I :nu0 :nuQ :nuref 
			:eta :K1 :K2
			:ntheta :nphi :Inu0 :nnu
			:intensity :expo
;		smoothedIntensity=:intensity
		smoothedIntensity=GaussianSmooth 
			:nnu :smoothwidth :compwidth :intensity
		mi=(max :smoothedIntensity)
		if :mi != 0
		[	intensityScaled=:smoothedIntensity/:mi*I0*400-maxy
		]
		hideTurtle
		clearScreen
		WindowMode
		setpensize 0
		penup
		setXY -400 :intensityScaled.1
		pendown
		setXY :x :intensityScaled

		specxsc=(:specx+:xmax)/w*maxx*2-maxx
		setPC HSB 120 1 1
		penup
		setXY :specxsc.1 :specy.1
		pendown
		setXY :specxsc :specy
		penup
		setPC 0

		specxsc2=(:specx2+:xmax)/w*maxx*2-maxx
		setPC HSB 0 1 1
		penup
		setXY :specxsc2.1 :specy2.1
		pendown
		setXY :specxsc2 :specy2
		penup
		setPC 0

		markerxsc=(:markerx+:xmax)/w*maxx*2-maxx
		setPC HSB 240 1 1
		penup
		setXY :markerxsc -:maxy
		pendown
		setXY :markerxsc :maxy
		penup
		setPC 0

		setY -:maxy
		for [x -:ticksx 2*:ticksx]
		[	hx=:wq*(:x/:ticksx-0.5)
			hx=:hx+(remainder (:xmax+:xmin)/2 :wq/:ticksx)
			rx=:hx-(:xmax+:xmin)/2
			hx=:hx*:maxx*2/:w
			if :hx > :maxx or2 :hx < -:maxx [continueLoop]
			setX :hx
			PenDown
			setY -:maxy+:tick
			setY -:maxy
			PenUp
			setY :maxy-:tick
			PenDown
			setY :maxy
			PenUp
			setY -:maxy-:txtx
			seth 90
			label :rx
			setY -:maxy
		]

		setXY 0 -250
		Label "|B (T)|

		setXY 0 -280
		Label :p

		updateGraph
		if Key? 
		[	ch=readChar
			if :ch==Char 27 [break]
			if :ch=="s 
			[	penup
				setpensize 0
			;	boundingbox
				catch "error
				[	saveScreenVector "nmr_carterspectrum.pdf
				]
				err=error
				if :err != []
				[	print :err.2
				]
			]
		]
		(GC "true)
	]
end

to computeIntensityRatios
	r=(FloatArray 2*:I -:I+1/2)
	j=1
	ratios={{1}@0 {4 3}@0 {9 8 5}@0 {16 15 12 7}@0 {25 24 12 16 9}@0}
	for [m -:I+1 :I]
	[	j=Int :m-1/2
		factor=:ratios.(Int :I+1/2).abs :j
		intensity=FloatArray :ni
		intensity=NMRCarterIntensity_m
			:I :nu0 :nuQ :nuref 
			:eta :K1 :K2
			:ntheta :nphi :Inu0 :nnu
			:intensity :m
		smoothedIntensity=GaussianSmooth :nnu :smoothwidth :intensity
		mi=(max :smoothedIntensity)
		if :mi != 0
		[	intensityScaled=:smoothedIntensity/:mi*factor*40-200
		]
		penup
		setpc HSB 180+180*:m/:I 1 1
		setXY -400 :intensityScaled.1
		pendown
		setXY :x :intensityScaled
		setitem :j :r 0+:intensity
	]
;		show :r
end

to readfile f
	local [x y d]
	x=[]
	y=[]
	openRead :f
	setReader :f
	while [not eof?]
	[	d=readlist
		x=fput :d.1 :x
		y=fput :d.2 :y
	]
	setReader []
	close :f
	output list reverse FloatArray :x reverse FloatArray :y
end