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