#pragma rtGlobals=3		// Use modern global access method and strict wave access.

Function init_angle(angle_res)
	Variable angle_res	
	Variable i
	Make/N=(360/angle_res)/D angle
	for (i=1;i<(numpnts(angle));i+=1)
		angle[i]=angle[i-1]+angle_res
	EndFor
	Make/N=(8)/D angle_ticklocator
	for (i=1;i<(numpnts(angle_ticklocator));i+=1)
		angle_ticklocator[i]=angle_ticklocator[i-1]+45
	EndFor
	angle_ticklocator=angle_ticklocator/angle_res
	Make/O/T angle_ticklabel = {"N", "NE", "E", "SE", "S", "SW", "W", "NW"}
End Function

Function init_speed(max_speed,speed_res)
	Variable max_speed,speed_res	
	Variable i
	Make/N=(max_speed/speed_res)/D speed
	for (i=1;i<(numpnts(speed));i+=1)
		speed[i]=speed[i-1]+speed_res
	EndFor
End Function

Function init_mx()
	wave angle,speed
	Make/N=(numpnts(angle),numpnts(speed))/D mx
	mx=nan
	//wave angle_ticklocator,angle_ticklabel
	//Display;AppendImage mx
	//ModifyGraph userticks(bottom)={angle_ticklocator,angle_ticklabel}
	//ModifyImage mx ctab= {*,*,Rainbow,1}
	//ModifyGraph axisEnab(bottom)={0,0.9}
	//ColorScale/C/N=text0/F=0/A=MC  ctab={0,100,Rainbow,1}
	//ColorScale/C/N=text0/A=RC/X=-1.65/Y=2.84
End Function

Function Kernel1(x)
	variable x
	return exp(-0.5*x*x)/sqrt(2*pi)
End Function

Function Kernel2(x)
	variable x
	if (0.75*(1-(x*x))<0)
		return 0
	Else
		return 0.75*(1-(x*x))
	EndIf
End Function

Function NWR_init()
	wave ws, wd
	
	duplicate ws diff_angle diff_speed angle_Kernel1 speed_Kernel2, Ktot
	diff_angle=nan
	diff_speed=nan
	angle_Kernel1=nan
	speed_Kernel2=nan
	Ktot=nan
End Function

Function angle_vec_diff(a1,a2)
	variable a1,a2
	variable b1,b2
	b1=a1*pi/180
	b2=a2*pi/180
	variable c
	c=cos(b1)*cos(b2)+sin(b1)*sin(b2)
	return acos(c)*180/pi
End Function

Function NWR(sigma,h)
	variable sigma, h
	wave ws, wd, angle, speed, conc, mx, diff_angle, diff_speed, angle_Kernel1, speed_Kernel2, Ktot
	variable i
	for (i=0;i<(numpnts(angle));i+=1)
		variable j
		for (j=0;j<(numpnts(speed));j+=1)
			diff_angle=angle_vec_diff(angle[i],wd)/sigma
			angle_Kernel1=Kernel1(diff_angle)
			diff_speed=(speed[j]-ws)/h
			speed_Kernel2=Kernel2(diff_speed)
			Ktot=speed_Kernel2*angle_Kernel1
			MatrixTranspose Ktot
			MatrixMultiply Ktot,conc
			wave M_product
			mx[i][j]=M_product[0][0]/sum(Ktot)
		EndFor		
	EndFor
End Function