This routine produces topographic map style plots from a 2D array representing
measured or calculated values on a regular rectangular grid.

Contours are obtained by inserting contrasting colours in the color table at
a specified interval.

Although the resulting contours are sometimes not perfect lines, this simple
technique may be useful for many practical purposes.

The routine was written for PBCC.

[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]'                     SHADED / CONTOUR PLOTS FROM 2D ARRAYS[/FONT]
[FONT=Courier New]'                                   for PBCC[/FONT]
[FONT=Courier New]'                          Arie Verheul - october 2008[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]' This routine produces topographic map style plots from a 2D array representing[/FONT]
[FONT=Courier New]' measured or calculated values on a regular rectangular grid.[/FONT]
[FONT=Courier New]' Contours are obtained by inserting contrasting colours in the color table at[/FONT]
[FONT=Courier New]' a specified interval.[/FONT]
[FONT=Courier New]' Although the resulting contours are sometimes not perfect lines, this simple[/FONT]
[FONT=Courier New]' technique may be useful for many practical purposes.[/FONT]
[FONT=Courier New]' The interpolated array ScreenArray() may be used to retrieve local values.[/FONT]
[FONT=Courier New]' For clarity this option, and some control options, were omitted from the code.[/FONT]
[FONT=Courier New]' The color range may be limited according to a specific interest.[/FONT]
[FONT=Courier New]' Therefore the color scale is extended at both sides with a dark green overflow color.[/FONT]
[FONT=Courier New]' The "zero" contour is dark gray.[/FONT]
[FONT=Courier New]' To use the plot routine in a different context requires the Subs DetermineExtremes,[/FONT]
[FONT=Courier New]' MapArrayToWindow and ColorPlot, plus the Type variable ContourParmset.[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]' For those interested in details about the example used here:[/FONT]
[FONT=Courier New]' Finite difference is a simple algorithm to (approximately) solve the Laplacian.[/FONT]
[FONT=Courier New]' It uses a 2D grid of points, usually representing some aspect of a physical structure.[/FONT]
[FONT=Courier New]' If grid spacing in both directions is equal, the algorithm involves repeatedly setting[/FONT]
[FONT=Courier New]' each gridpoint (except for boundary points) at the average of its four neighbours.[/FONT]
[FONT=Courier New]' The example uses a grid of 101 x 101 gridpoints.[/FONT]
[FONT=Courier New]' The plot shows an example of "field-mapping", a technique to use the field[/FONT]
[FONT=Courier New]' from point or line sources (which can be easily calcuated) to model something else.[/FONT]
[FONT=Courier New]' The idea behind this is that each equipotential contour may be interpreted[/FONT]
[FONT=Courier New]' as a conducting surface (or an electrode).[/FONT]
[FONT=Courier New]' Therefore, the calculated 2D field from five line sources perpendicular to the[/FONT]
[FONT=Courier New]' plane may as well represent the field from five rods.[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]#Compiler PBCC 4.04[/FONT]
[FONT=Courier New]#Compile Exe[/FONT]
[FONT=Courier New]#Dim All[/FONT]
[FONT=Courier New]#Console Off[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Type ContourParmset[/FONT]
[FONT=Courier New]  NumColor  As Long       ' Number of shades in positive and negative half of range[/FONT]
[FONT=Courier New]  PlotRange As Single     ' Maximum value to be rendered (positive value)[/FONT]
[FONT=Courier New]  ContInt   As Single     ' Contour interval, in same units as PlotRange (positive value)[/FONT]
[FONT=Courier New]  ArrayMax  As Single     ' Maximum and minimum values, determined internally[/FONT]
[FONT=Courier New]  ArrayMin  As Single[/FONT]
[FONT=Courier New]End Type[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Declare Sub      FiniteDifference   (PlotArray()   As Single)[/FONT]
[FONT=Courier New]Declare Sub      DetermineExtremes  (PlotArray()   As Single, PlotParms As ContourParmset)[/FONT]
[FONT=Courier New]Declare Sub      MapArrayToWindow   (PlotArray()   As Single, ScreenArray() As Single)[/FONT]
[FONT=Courier New]Declare Sub      ColorPlot          (ScreenArray() As Single, PlotParms As ContourParmset)[/FONT]
[FONT=Courier New]Declare Function IsWindow Lib "USER32.DLL" Alias "IsWindow" (ByVal hWnd As Dword) As Long[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Function PBMain () As Long[/FONT]
[FONT=Courier New]  Local hWnd As Dword[/FONT]
[FONT=Courier New]  Dim ScreenArray ()  As Single[/FONT]
[FONT=Courier New]  Dim PlotArray()     As Single[/FONT]
[FONT=Courier New]  Dim PlotParms       As ContourParmset[/FONT]
[FONT=Courier New]  ' Setup window[/FONT]
[FONT=Courier New]  Graphic Window "Electric potential plot",200,100,600,600 To hWnd[/FONT]
[FONT=Courier New]  Graphic Attach hWnd, 0[/FONT]
[FONT=Courier New]  Graphic Color  %White, %Black[/FONT]
[FONT=Courier New]  Graphic Clear[/FONT]
[FONT=Courier New]  ' Plot settings (may be changed)[/FONT]
[FONT=Courier New]  PlotParms.NumColor  = 250[/FONT]
[FONT=Courier New]  PlotParms.PlotRange = 500[/FONT]
[FONT=Courier New]  PlotParms.ContInt   =  50[/FONT]
[FONT=Courier New]  ' Main program[/FONT]
[FONT=Courier New]  FiniteDifference PlotArray()                ' Produce some data to be plotted[/FONT]
[FONT=Courier New]  DetermineExtremes PlotArray(), PlotParms    ' Determine Max/Min values[/FONT]
[FONT=Courier New]  MapArrayToWindow PlotArray(), ScreenArray() ' Calculate interpolated values for each pixel[/FONT]
[FONT=Courier New]  ColorPlot ScreenArray(), PlotParms          ' Translate values to shades[/FONT]
[FONT=Courier New]  While IsWindow(hWnd): Sleep 1000: Wend      ' Terminate program after window closure[/FONT]
[FONT=Courier New]End Function[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Sub FiniteDifference (PlotArray() As Single)[/FONT]
[FONT=Courier New]  ' The purpose of this Sub is just to produce some data to be plotted[/FONT]
[FONT=Courier New]  ' It must not be considered as a good or general example to use this algorithm[/FONT]
[FONT=Courier New]  %NumInt = 100[/FONT]
[FONT=Courier New]  Local N, M, I As Long[/FONT]
[FONT=Courier New]  Local XS1, XS2, XS3, YS1, YS2, YS3 As Long[/FONT]
[FONT=Courier New]  ReDim PlotArray(%NumInt,%NumInt) As Single[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' Set internal boundaries[/FONT]
[FONT=Courier New]  XS1 = .2 * %NumInt              ' Location of internal boundary points[/FONT]
[FONT=Courier New]  XS2 = %NumInt - XS1[/FONT]
[FONT=Courier New]  YS1 = .2 * %NumInt[/FONT]
[FONT=Courier New]  YS2 = %NumInt - YS1[/FONT]
[FONT=Courier New]  XS3 = %NumInt / 2[/FONT]
[FONT=Courier New]  YS3 = %NumInt / 2[/FONT]
[FONT=Courier New]  PlotArray(XS1,YS1) =   1250     ' Values assigned to internal boundary points[/FONT]
[FONT=Courier New]  PlotArray(XS2,YS2) =   1250[/FONT]
[FONT=Courier New]  PlotArray(XS1,YS2) = - 1375[/FONT]
[FONT=Courier New]  PlotArray(XS2,YS1) = - 1375[/FONT]
[FONT=Courier New]  PlotArray(XS3,YS3) =   1125[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  For I = 1 To 2000                  ' Iterations[/FONT]
[FONT=Courier New]      For M = 1 To %NumInt - 1       ' Vert. direction[/FONT]
[FONT=Courier New]          For N = 1 To %NumInt - 1   ' Hor. direction[/FONT]
[FONT=Courier New]              If M = YS1 Or M = YS2 Then[/FONT]
[FONT=Courier New]                  If N = XS1 Or N = XS2 Then Iterate For  ' Skip internal boundary points[/FONT]
[FONT=Courier New]              End If[/FONT]
[FONT=Courier New]              If M = YS3 And N = XS3 Then Iterate For[/FONT]
[FONT=Courier New]              PlotArray(N,M) = (PlotArray(N-1,M) + PlotArray(N+1,M)_[/FONT]
[FONT=Courier New]                               + PlotArray(N,M-1) + PlotArray(N,M+1)) / 4[/FONT]
[FONT=Courier New]          Next[/FONT]
[FONT=Courier New]      Next[/FONT]
[FONT=Courier New]  Next[/FONT]
[FONT=Courier New]End Sub[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Sub DetermineExtremes (PlotArray() As Single, PlotParms As ContourParmset)[/FONT]
[FONT=Courier New]  ' Determines Maximum and Minimum value in Plotarray()[/FONT]
[FONT=Courier New]  Local M, N As Long[/FONT]
[FONT=Courier New]  Local TestValue, ArrayMax, ArrayMin As Single[/FONT]
[FONT=Courier New]  ArrayMin = PlotArray(0,0)[/FONT]
[FONT=Courier New]  ArrayMax = ArrayMin[/FONT]
[FONT=Courier New]  For M = 0 To UBound(PlotArray, 2)[/FONT]
[FONT=Courier New]      For N = 0 To UBound(PlotArray, 1)[/FONT]
[FONT=Courier New]          TestValue = PlotArray(N,M)[/FONT]
[FONT=Courier New]          'Skip poles, if present[/FONT]
[FONT=Courier New]          If (TestValue > 1E35) Or (TestValue < -1E35) Then Iterate For[/FONT]
[FONT=Courier New]          ArrayMax = Max(TestValue, ArrayMax)[/FONT]
[FONT=Courier New]          ArrayMin = Min(TestValue, ArrayMin)[/FONT]
[FONT=Courier New]      Next[/FONT]
[FONT=Courier New]  Next[/FONT]
[FONT=Courier New]  PlotParms.ArrayMax = ArrayMax[/FONT]
[FONT=Courier New]  PlotParms.ArrayMin = ArrayMin[/FONT]
[FONT=Courier New]End Sub[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Sub MapArrayToWindow (PlotArray() As Single, ScreenArray() As Single)[/FONT]
[FONT=Courier New]  ' This sub calculates linearly interpolated values from PlotArray()[/FONT]
[FONT=Courier New]  ' for each pixel in ScreenArray()[/FONT]
[FONT=Courier New]  ' Note: in the following the rectangular area enclosed by[/FONT]
[FONT=Courier New]  ' four neighbouring gridpoints is designated as an element[/FONT]
[FONT=Courier New]  Local WindowWidth, WindowHeight As Long[/FONT]
[FONT=Courier New]  Local Xmax, Ymax As Long[/FONT]
[FONT=Courier New]  Local I, ElementRow, Pix, PixRow, CurrentRow As Long[/FONT]
[FONT=Courier New]  Local X, PixRowPos, VertIntFactor As Single[/FONT]
[FONT=Courier New]  Graphic Get Client To WindowWidth, WindowHeight[/FONT]
[FONT=Courier New]  Xmax = UBound(PlotArray(),1)[/FONT]
[FONT=Courier New]  Ymax = UBound(PlotArray(),2)[/FONT]
[FONT=Courier New]  Dim TopRow(Xmax)               As Single   ' Top nodes of current element row[/FONT]
[FONT=Courier New]  Dim BotRow(Xmax)               As Single   ' Bottom nodes of current element row[/FONT]
[FONT=Courier New]  Dim RowDiff(Xmax)              As Single   ' Differential (TopRow() - BotRow)[/FONT]
[FONT=Courier New]  Dim InterpRow(Xmax)            As Single   ' Interpolated values for position of pixel row[/FONT]
[FONT=Courier New]  Dim ElDiff(Xmax - 1)           As Single   ' Hor. differential in InterpRow()[/FONT]
[FONT=Courier New]  Dim HorIntNumber (WindowWidth) As Long     ' Element to use for interpolation for specific pixel[/FONT]
[FONT=Courier New]  Dim HorIntFactor (WindowWidth) As Single   ' Interpolation factor for specific pixel[/FONT]
[FONT=Courier New]  ReDim ScreenArray (1 To WindowWidth, 1 To WindowHeight) As Single[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' Make tables for pixel locations and interpolation factors in element grid coordinates[/FONT]
[FONT=Courier New]  ' in horizontal direction; these tables are valid for each row of pixels.[/FONT]
[FONT=Courier New]  For Pix = 1 To WindowWidth          ' For each pixel in horizontal row[/FONT]
[FONT=Courier New]      X = Pix * Xmax / WindowWidth    ' > Pixel position in element grid coordinates[/FONT]
[FONT=Courier New]      HorIntNumber (Pix) = Int  (X)   ' > Element where pixel is located[/FONT]
[FONT=Courier New]      HorIntFactor (Pix) = Frac (X)   ' > Interpolation factor within element[/FONT]
[FONT=Courier New]  Next[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' Do the actual interpolation, working row by row[/FONT]
[FONT=Courier New]  CurrentRow = -1                                  ' Counter to keep track of current element row[/FONT]
[FONT=Courier New]  For PixRow = 1 To WindowHeight                   ' Work row by row[/FONT]
[FONT=Courier New]      PixRowPos = Ymax * PixRow / WindowHeight     ' Vertical position of pixel row[/FONT]
[FONT=Courier New]                                                   ' in element grid coordinates[/FONT]
[FONT=Courier New]      ElementRow    = Int  (PixRowPos)             ' Element row to use for interpolation[/FONT]
[FONT=Courier New]      VertIntFactor = Frac (PixRowPos)             ' Interpolation factor[/FONT]
[FONT=Courier New]      If ElementRow <> CurrentRow Then             ' If interpolation has reached next[/FONT]
[FONT=Courier New]                                                   ' element row, data must be refreshed[/FONT]
[FONT=Courier New]          ' Full row of data is extracted from PlotArray() by placing 1D array on top of it[/FONT]
[FONT=Courier New]          ' Using MAT these data are manipulated with a full row at the time[/FONT]
[FONT=Courier New]          ReDim TopRow(Xmax) As Single At VarPtr (PlotArray(0, ElementRow + 1))[/FONT]
[FONT=Courier New]          ReDim BotRow(Xmax) As Single At VarPtr (PlotArray(0, ElementRow))[/FONT]
[FONT=Courier New]          Mat   RowDiff() = TopRow() - BotRow()    ' Vertical differentials[/FONT]
[FONT=Courier New]          CurrentRow = ElementRow[/FONT]
[FONT=Courier New]      End If[/FONT]
[FONT=Courier New]      Mat    InterpRow () = (VertIntFactor) * RowDiff()[/FONT]
[FONT=Courier New]      Mat    InterpRow () = BotRow() + InterpRow()[/FONT]
[FONT=Courier New]      ' InterpRow() now contains an interpolated row of data from PlotArray()[/FONT]
[FONT=Courier New]      ' at PixRowPos, the position of the row of pixels that is treated currently.[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]      ' Next within this row the differentials are calculated, and stored in ElDiff()[/FONT]
[FONT=Courier New]      For I = 0 To Xmax - 1           ' Horizontal differential in InterpRow()[/FONT]
[FONT=Courier New]          ElDiff(I) = InterpRow(I + 1) - InterpRow(I)[/FONT]
[FONT=Courier New]      Next[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]      ' From InterpRow() and ElDiff() the values for each pixel are calculated[/FONT]
[FONT=Courier New]      For Pix = 1 To WindowWidth      ' Work along a row of pixels[/FONT]
[FONT=Courier New]          I = HorIntNumber(Pix)       ' Element where pixel is located[/FONT]
[FONT=Courier New]          ScreenArray(Pix, PixRow) = InterpRow(I) + ElDiff(I) * HorIntFactor(Pix)[/FONT]
[FONT=Courier New]      Next        ' Next pixel[/FONT]
[FONT=Courier New]  Next            ' Next row of pixels[/FONT]
[FONT=Courier New]End Sub[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]Sub ColorPlot (ScreenArray() As Single, PlotParms As ContourParmset)[/FONT]
[FONT=Courier New]  ' This sub translates values associated with pixels into shades[/FONT]
[FONT=Courier New]  Local R,G,B As Integer[/FONT]
[FONT=Courier New]  Local N, M, ColorIndex As Long[/FONT]
[FONT=Courier New]  Local ColorScaleUBound, ColorScaleLBound As Long[/FONT]
[FONT=Courier New]  Local WindowWidth, WindowHeight As Long[/FONT]
[FONT=Courier New]  Local X, NormalisedIndex, NormalisedIndexSqr As Single[/FONT]
[FONT=Courier New]  Local WindowMap As String[/FONT]
[FONT=Courier New]  Dim   PixelPtr As Long Pointer[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' Set default values if necessary[/FONT]
[FONT=Courier New]  If PlotParms.NumColor <= 0 Then PlotParms.NumColor = 300[/FONT]
[FONT=Courier New]  PlotParms.PlotRange    = Abs(PlotParms.PlotRange)[/FONT]
[FONT=Courier New]  If PlotParms.PlotRange = 0 Then[/FONT]
[FONT=Courier New]     PlotParms.PlotRange = Max(Abs(PlotParms.ArrayMax), Abs(PlotParms.ArrayMin))[/FONT]
[FONT=Courier New]  End If[/FONT]
[FONT=Courier New]  PlotParms.ContInt      = Abs(PlotParms.ContInt)[/FONT]
[FONT=Courier New]  If PlotParms.ContInt   = 0 Then PlotParms.ContInt = PlotParms.PlotRange / 20[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ColorScaleUBound =  PlotParms.NumColor + 1[/FONT]
[FONT=Courier New]  ColorScaleLBound = -PlotParms.NumColor - 1[/FONT]
[FONT=Courier New]  Dim Shade(ColorScaleLBound To ColorScaleUBound) As Long[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' make color table, involving the following steps:[/FONT]
[FONT=Courier New]  ' a. a fractional value <NormalisedIndex> is calculated from the color index.[/FONT]
[FONT=Courier New]  ' b. using the expression contained in ColorAlgorithm the "positive" shading.[/FONT]
[FONT=Courier New]  '    colors are calculated (BGR format).[/FONT]
[FONT=Courier New]  ' c. swapping R and B value gives "negative" shading colors.[/FONT]
[FONT=Courier New]  ' d. contour colors are obtained by adding .35 to <NormalisedIndex>.[/FONT]
[FONT=Courier New]  '    which shifts them to a brighter value.[/FONT]
[FONT=Courier New]  ' e. if <NormalisedIndex> > 1 then white is used for contours.[/FONT]
[FONT=Courier New]  ' f. "zero" contour is assigned dark gray.[/FONT]
[FONT=Courier New]  ' g. at both sides of the color scale an overflow color, dark green, is added.[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' background colors[/FONT]
[FONT=Courier New]  For N = 1 To PlotParms.NumColor[/FONT]
[FONT=Courier New]      NormalisedIndex = N / PlotParms.NumColor[/FONT]
[FONT=Courier New]      GoSub ColorAlgorithm[/FONT]
[FONT=Courier New]  Next[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' contour colors[/FONT]
[FONT=Courier New]  For X = 0 To PlotParms.PlotRange Step PlotParms.ContInt[/FONT]
[FONT=Courier New]      NormalisedIndex = X / PlotParms.PlotRange[/FONT]
[FONT=Courier New]      N               = NormalisedIndex * PlotParms.NumColor[/FONT]
[FONT=Courier New]      NormalisedIndex = NormalisedIndex + .35[/FONT]
[FONT=Courier New]      If NormalisedIndex <= 1 Then    ' the above line may produce values > 1[/FONT]
[FONT=Courier New]          GoSub ColorAlgorithm[/FONT]
[FONT=Courier New]      Else[/FONT]
[FONT=Courier New]          Shade (N) = %White          ' use white instead[/FONT]
[FONT=Courier New]          Shade(-N) = %White[/FONT]
[FONT=Courier New]      End If[/FONT]
[FONT=Courier New]  Next[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  Shade(0) = RGB(80,80,80)                      ' "zero" contour[/FONT]
[FONT=Courier New]  Shade( PlotParms.NumColor + 1) = RGB(0,80,0)  ' overflow color[/FONT]
[FONT=Courier New]  Shade(-PlotParms.NumColor - 1) = RGB(0,80,0)[/FONT]
[FONT=Courier New]  ' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]  ' GRAPHIC GET BITS / GRAPHIC SET BITS is used here because it is[/FONT]
[FONT=Courier New]  ' much faster than GRAPHIC SET PIXEL[/FONT]
[FONT=Courier New]  Graphic Get Client To WindowWidth, WindowHeight[/FONT]
[FONT=Courier New]  Graphic Get Bits To WindowMap       ' Get copy of video memory[/FONT]
[FONT=Courier New]  PixelPtr = StrPtr(WindowMap) + 8    ' First two Long Int. contain size[/FONT]
[FONT=Courier New]  For M = WindowHeight To 1 Step -1   ' Origin at left bottom corner[/FONT]
[FONT=Courier New]                                      ' For origin at left top corner reverse scanning direction[/FONT]
[FONT=Courier New]      For N = 1 To WindowWidth[/FONT]
[FONT=Courier New]      ColorIndex = ScreenArray(N,M) * (PlotParms.NumColor/PlotParms.PlotRange)[/FONT]
[FONT=Courier New]      ' overflow check[/FONT]
[FONT=Courier New]      ColorIndex = Min(ColorIndex, ColorScaleUBound)[/FONT]
[FONT=Courier New]      ColorIndex = Max(ColorIndex, ColorScaleLBound)[/FONT]
[FONT=Courier New]      @PixelPtr = Shade(ColorIndex)   ' Write color to Windowmap[/FONT]
[FONT=Courier New]      Incr PixelPtr[/FONT]
[FONT=Courier New]      Next[/FONT]
[FONT=Courier New]  Next[/FONT]
[FONT=Courier New]  Graphic Set Bits WindowMap          ' Put modified WindowMap back[/FONT]
[FONT=Courier New]Exit Sub[/FONT]
[FONT=Courier New]' ---------------------------------------------------------------------------[/FONT]
[FONT=Courier New]ColorAlgorithm:[/FONT]
[FONT=Courier New]  NormalisedIndexSqr = NormalisedIndex^2[/FONT]
[FONT=Courier New]  R = 600 * NormalisedIndex - 345 * NormalisedIndexSqr   ' color expression[/FONT]
[FONT=Courier New]  G = 255 * NormalisedIndexSqr[/FONT]
[FONT=Courier New]  B =   0[/FONT]
[FONT=Courier New]  Shade (N) = RGB(B,G,R)              ' BGR color format[/FONT]
[FONT=Courier New]  Shade(-N) = RGB(R,G,B)              ' for "negative" colors swap R and B[/FONT]
[FONT=Courier New]Return[/FONT]
[FONT=Courier New]End Sub[/FONT]