; requires co-ordinates of the centre in data terms, with datum
; SpecUnitArea is an option: set to the number of years (x-axis) representing the area when Datum=1
; returns the vectors describing a circle in data terms

pro CircleRefs, CentreX,CentreY, AxisLenX,AxisLenY, Datum, VecX,VecY, SpecUnitArea=SpecUnitArea

VecX = dblarr (361)
VecY = dblarr (361)

if keyword_set(SpecUnitArea) then begin
  UnitArea = double(SpecUnitArea)
endif else begin
  UnitArea = double(AxisLenX) / 1000.0
endelse

Area   = double(Datum) * UnitArea
Radius = sqrt (Area/3.1415927)
Ratio  = double(AxisLenY) / double(AxisLenX)

for XAngle = 0, 89 do begin
    Opposite = double(Radius * sin ((89.5-double(XAngle)) * 3.1415927 / 180.0))
    Adjacent = double(sqrt (Radius*Radius - Opposite*Opposite))
    
    VecX (  0+XAngle) = CentreX + Adjacent
    VecY (  0+XAngle) = CentreY + Opposite * Ratio
endfor

for XAngle = 0, 89 do begin
    Opposite = double(Radius * sin ((double(XAngle)+0.5) * 3.1415927 / 180.0))
    Adjacent = double(sqrt (Radius*Radius - Opposite*Opposite))
    
    VecX ( 90+XAngle) = CentreX + Adjacent
    VecY ( 90+XAngle) = CentreY - Opposite * Ratio
endfor

for XAngle = 0, 89  do begin
    Opposite = double(Radius * sin ((89.5-double(XAngle)) * 3.1415927 / 180.0))
    Adjacent = double(sqrt (Radius*Radius - Opposite*Opposite))
    
    VecX (180+XAngle) = CentreX - Adjacent
    VecY (180+XAngle) = CentreY - Opposite * Ratio
endfor

for XAngle = 0, 89 do begin
    Opposite = double(Radius * sin ((double(XAngle)+0.5) * 3.1415927 / 180.0))
    Adjacent = double(sqrt (Radius*Radius - Opposite*Opposite))
    
    VecX (270+XAngle) = CentreX - Adjacent
    VecY (270+XAngle) = CentreY + Opposite * Ratio
endfor

VecX(360) = VecX(0) & VecY(360) = VecY(0)

end