Taylor diagrams

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

Macro RedrawTaylorDiagramTicksGrids()
	RedrawCorrelationTicks()
	RedrawRadiusGrids()
	RedrawIsoLines()
End

Proc RedrawCorrelationTicks()
	Variable labelFontSize=12
	Variable drawMajorCorrelationsGrid=1
	Variable majorCorrelationsTo=0.25
	TopTaylorDiagram(1)
	GetAxis/Q bottom
	Variable ticksAtRadius= V_Max

	fReDrawCorrelationTicks(ticksAtRadius, labelFontSize, drawMajorCorrelationsGrid, majorCorrelationsTo)
End

Proc RedrawRadiusGrids()
	Variable firstRadius = 0.25
	Variable radiusInc = 0.25
	Variable maxRadius = 1.5
	TopTaylorDiagram(1)
	fRedrawRadiusGrids(firstRadius, radiusInc, maxRadius)
End

Proc RedrawIsoLines()
	Variable centerX= 1.15
	Variable centerY= 0
	Variable isoRadiusIncrement= 0.4
	Variable numberOfLines= 4
	TopTaylorDiagram(1)
	fRedrawIsoLines(centerX, centerY, isoRadiusIncrement,numberOfLines)
End

Macro ClearTaylorDiagramTicksAndGrids()

	TopTaylorDiagram(1)
	SetDrawEnv push
	SetDrawLayer/K ProgAxes	// IsoLines
	SetDrawLayer/K ProgBack	// RadiusGrids
	SetDrawLayer/K UserAxes	// CorrelationTicks
	SetDrawEnv pop
End

Function/S TopTaylorDiagram(bringToFront)
	Variable bringToFront
	
	String taylorGraphs= WinList("TaylorDiagram*",";","WIN:1,VISIBLE:1")
	String graphName=StringFromList(0,taylorGraphs)	// first one is the topmost one
	if( strlen(graphName) && bringToFront )
		DoWindow/F $graphName
	endif
	return graphName
End


Function fRedrawIsoLines(centerX, centerY, isoRadiusIncrement,numberOfLines)
	Variable centerX, centerY, isoRadiusIncrement,numberOfLines

	SetDrawEnv push
	SetDrawLayer/K ProgAxes	// Note: different layer than fReDrawCorrelationTicks and fRedrawRadiusGrids
	SetDrawEnv xcoord= bottom,ycoord= left,dash= 0,linefgc=(0,65535,0), fillpat=0,save	// solid, green
	
	// indicate center with a drawn marker
	SetDrawEnv textxjust=1, textyjust=1, textrgb=(0,65535,0)
	DrawText centerX, centerY,"o"
	GetAxis/Q bottom
	Variable radiusMax= V_Max
	Variable i
	for(i=0; i<numberOfLines; i+=1)
		Variable radius= (i+1)*isoRadiusIncrement
		Variable startDegrees=StartAngleFor(centerX, radius, radiusMax)
		Variable endDegrees=StopAngleFor(centerX, radius, radiusMax)
		DrawArc/X/Y centerX, centerY, radius, startDegrees, endDegrees

		Variable labelX, labelY
		FindIntersectionOfArcWithRadius(centerX,radius, labelX, labelY)
		Variable labelRadians= atan2(labelY,centerx-labelx)
		Variable labelDegrees= 90-labelRadians * 180 / pi
		SetDrawEnv textrot=labelDegrees, textxjust=1,textyjust=1,textrgb=(0,65535,0)
		String labelText= num2str(radius)
		DrawText labelX, labelY,labelText
	endfor
	
	SetDrawEnv pop
End

// put label at the intersection of the iso circle
// and a circle centered at the origin with radiusFromZero
Function FindIntersectionOfArcWithRadius(radiusFromZero,radiusFromIso, labelX, labelY)
	Variable radiusFromZero, radiusFromIso
	Variable &labelX, &labelY	// outputs
	
	// math.stackexchange.com/questions/256100/how-can-i-find-the-points-at-which-two-circles-intersect
	Variable x1= 0, y1=0, r1=radiusFromZero				// circle1 centered at the origin with radiusFromZero
	Variable x2= radiusFromZero, y2=0, r2= radiusFromIso	// circle2 at radiusFromZero,0 with radiusFromIso
	
	// -2x(x1-x2) - 2y(y1-y2) = (r1^2-r2^2) - (x1^2 - x2^2) - (y1^2 - y2^2)
	// in our case both y1 and y2 are zero, so we solve for x:
	// -2x(x1-x2) = (r1^2-r2^2) - (x1^2 - x2^2)
	//
	// or x = [(r1^2-r2^2) - (x1^2 - x2^2)]/(2x2-2x1)
	labelX= ((r1*r1-r2*r2) - (x1*x1 - x2*x2))/(2*x2-2*x1)
	labelY= sqrt(radiusFromZero*radiusFromZero - labelX*labelX)
End
	

// returns degrees, presumes the graph's angle range is 0-90 degrees, and the arc center is on the x axis between 0 and radiusMax
Function StartAngleFor(centerX, radius, radiusMax)
	Variable centerX, radius, radiusMax
	
	Variable angleInDegrees= 0
	if( (centerX + radius) > radiusMax )
		// intersects along radiusMax, use Law of Cosines since we know all three side lengths
		Variable numerator= centerX*centerX+radius*radius-radiusMax*radiusMax
		Variable denominator= 2*centerX*radius
		Variable ratio= numerator/denominator
		Variable complementInRadians= acos(ratio)
		angleInDegrees= 180/pi*(pi-complementInRadians)
	endif
	
	return angleInDegrees	
End

// returns degrees, presumes the graph's angle range is 0-90 degrees, and the arc center is on the x axis between 0 and radiusMax
Function StopAngleFor(centerX,  radius, radiusMax)
	Variable centerX,  radius, radiusMax
	
	Variable angleInDegrees= 180
	if( (centerX - radius) < 0 )
		// intersects along x=0
		Variable numerator= centerX
		Variable denominator= radius
		Variable ratio= numerator/denominator
		Variable complementInRadians= acos(ratio)
		angleInDegrees= 180/pi*(pi-complementInRadians)
	endif
	
	return angleInDegrees	
End


Function fRedrawRadiusGrids(firstRadius, radiusInc, maxRadius)
	Variable firstRadius, radiusInc, maxRadius
	
	SetDrawEnv push
	SetDrawLayer/K ProgBack	// Note: different layer than fReDrawCorrelationTicks and fRedrawIsoLines
	SetDrawEnv xcoord= bottom,ycoord= left,dash= 3,save

	// always draw a grid at radius=1 as solid black
	Variable radius= 1.0
	SetDrawEnv xcoord= bottom,ycoord= left,dash= 0,linefgc=(0,0,0),fillpat= 0, save
	DrawArc /X/Y 0,0,radius,0,90
	
	// Draw the other grids as dashed blue
	SetDrawEnv dash=3, linefgc=(32768,40777,65535),save	// light blue
	for( radius= firstRadius; radius <= maxRadius; radius += radiusInc )
		if( abs(radius - 1.0) > radiusInc/4 )	// don't draw over the black grid
			DrawArc /X/Y 0,0,radius,0,90
		endif
	endfor
	SetDrawEnv pop
End

Function fReDrawCorrelationTicks(ticksAtRadius, labelFontSize, drawMajorCorrelationsGrid, majorCorrelationsTo)
	Variable ticksAtRadius
	Variable labelFontSize
	Variable drawMajorCorrelationsGrid
	Variable majorCorrelationsTo
	
	// the following should perhaps be inputs instead of being hard-coded here
	Make/O majorTicks={0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99, 1.0 }
	Make/O minorTicks={0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.91, 0.92, 0.93, 0.94, 0.96, 0.97, 0.98}
	Variable majorTickLen = 0.04	// radius units
	Variable minorTickLen = majorTickLen * 0.5	// radius units

	SetDrawEnv push
	SetDrawLayer/K UserAxes
	SetDrawEnv xcoord= bottom,ycoord= left,dash= 0,fsize=labelFontSize,fillpat=0,save
	DrawArc /X/Y 0,0,ticksAtRadius,0,90	// draw axis

	Variable i, n, innerRadius, correlation, x1, x2, y1, y2, radiansForCorrelation, yScale, labelDegrees, labelRadius

	// draw major ticks
	n= numpnts(majorTicks)
	innerRadius= ticksAtRadius - majorTickLen
	for(i=0; i<n; i+=1 )
		correlation= majorTicks[i]
		x1 = innerRadius * correlation
		x2 = ticksAtRadius * correlation
		radiansForCorrelation = acos(correlation)
		yScale= sin(radiansForCorrelation)
		y1 = innerRadius * yScale
		y2 = ticksAtRadius * yScale
		if( drawMajorCorrelationsGrid && correlation != 0 && correlation != 1.0 )
			SetDrawEnv push
			SetDrawEnv dash=3, linefgc=(32768,40777,65535)	// light blue
			DrawLine majorCorrelationsTo*correlation,majorCorrelationsTo* yScale,x1,y1
			SetDrawEnv pop
		endif
		DrawLine x1,y1,x2,y2
		String labelText= num2str(correlation)
		labelDegrees= radiansForCorrelation * 180 / pi
		labelRadius= ticksAtRadius + majorTickLen
		SetDrawEnv push
		SetDrawEnv rotate=labelDegrees, textyjust=1
		DrawText labelRadius,0,labelText
		SetDrawEnv pop
	endfor

	// draw minor ticks
	n= numpnts(minorTicks)
	innerRadius= ticksAtRadius - minorTickLen
	for(i=0; i<n; i+=1 )
		correlation= minorTicks[i]
		x1 = innerRadius * correlation
		x2 = ticksAtRadius * correlation
		radiansForCorrelation = acos(correlation)
		yScale= sin(radiansForCorrelation)
		y1 = innerRadius * yScale
		y2 = ticksAtRadius * yScale
		DrawLine x1,y1,x2,y2
	endfor

	SetDrawEnv pop
End

Macro NewTaylorDiagram() : Graph
	PauseUpdate; Silent 1		// building window...
	Make/O/N=1 invisible=NaN
	Display /W=(239,135,706,576)/N=TaylorDiagram invisible as "Taylor Diagram"
	ModifyGraph margin(left)=64,margin(bottom)=47,margin(top)=41,margin(right)=50,width={Plan,1,bottom,left}
	ModifyGraph standoff=0
	ModifyGraph manTick(left)={0,0.25,0,2},manMinor(left)={0,50}
	ModifyGraph manTick(bottom)={0,0.25,0,2},manMinor(bottom)={0,50}
	SetAxis left 0,1.7
	SetAxis bottom 0,1.7
	TextBox/C/N=text0/O=-50/F=0/A=MC/X=32.29/Y=26.91 "Correlation"
EndMacro