' this is PhaseArray2D-01.bas   in FreeBasic code
' use the picture PhaseArray.gif in /PhaseArray/
'
'     the purpose of this version is to do 2D:
'         increase the resolution around beam center 0
'             hi scale to 0-2 degrees, low scale 2-90 degrees
'   enhancement over version 02 - after first pass with 
'       straight phasing, go to correlation to check noise floor
'
' game is to rotate a flat wave front through some max angle? 
'       centered at the antenna farm reference point
'
'  all dimensions are in wavelengths
'   plan a) set up operating parameters
'        b) swing arc, storing for plotting, or plotting as we go?
'           plot max y in graph delta x, expecting max to be # antennas,
'                    min to be ??? -1??
'  analysis and output
'    I think I like voltage better, but folks like power better 
'       for instance, the output of an FFT can be complex, but more
'          easily shown as power spectra - screwing up phase and reversibility :-|
'    And folks use beamwidth as what?
'              down 3 db voltage and half power??
'        so - beam width described below is 
'    Now, about logarithmic outputs?  y axis- 

' Lets do it in 0.5 degree increments, out to ??
'
' Well - interesting - in large arrays, ie ATA, the offaxis peaks
' might well be treated as an electrical noise??
'    Think I will enhance the program to treat the hard stuff to visualize
'       as noise, and report the average a ? db 
' 
'23456789012345678901234567890123456789012345678901234567890123456789012345678901

' DEFDBL A-Z   ' do all double, why muck around?  
DEFSNG A-Z   ' do all double, why muck around?  

CONST Pi = 3.1415926535897932384626433832795
CONST TwoPi = Pi*2
CONST Rad2Deg = 360/TwoPi
CONST Deg2Rad = TwoPi/360
CONST NumAnt = 350  ' Full ATA
CONST DiaWidthWave = 70000  ' width in wavelengths- expanded ATA?
'CONST NumAnt = 42  '  ATA phase 42
'CONST DiaWidthWave = 7000  ' width in wavelengths- phase 42 ATA?
DIM AntX(350)       '  location of antenna in X from ref, in units of wave lengths
DIM AntY(350)       '  location of antenna in Y from ref, in units of wave lengths
DIM AntPhase(350)  ' phase from wave plane - recalculated each different wave plane angle
DIM  PowerArrayByColor(10)
DIM MaskPowerArrayByColor(10)
For x = 1 to 10
    PowerArrayByColor(x) = 0
    MaskPowerArrayByColor(x) = 0
Next x


 SCREEN 21, 32, 1,0 ' 15 ' 14  ' 12 '13
    ' 16-512x384; 17-640x400; 18-640x480; 19-800x600; 20-1024x768; 21-1280x1024
         colorGreen = &H0000FF00
         colorRed   = &H00FF0000
         colorWhite = &H00FFFFFF
         colorBlue  = &H000000FF
         colorYellow= &H00FFFF00
         colorGray  = &H00AAAAAA
         colorViolet= &H0000A0A0
     
randomize(0)
n = 1
AntRadiusWave = DiaWidthWave/2

        AntX(1) = 0
        AntY(1) = 0
        AntX(2) = 0
        AntY(2) = AntRadiusWave
        AntX(3) = AntRadiusWave
        AntY(3) = 0
        AntX(4) = -AntRadiusWave
        AntY(4) = 0
        AntX(5) = 0
        AntY(5) = -AntRadiusWave

While n <= NumAnt 
    x = DiaWidthWave*rnd(1)-AntRadiusWave
    y = DiaWidthWave*rnd(1)-AntRadiusWave
    r = (x*x + y*y)^0.5
    If r < AntRadiusWave then ' in diameter of circle
        AntX(n) = x
        AntY(n) = y
        'Print "Ant";n, "X=";x, "Y=";y, "RadiusWave=";r
        n = n+1
    End If
Wend ' n
  WHILE INKEY$ = "": sleep (10):WEND

RadiusCirc = 300 ' pixels
RadiusField = 70 ' degrees

Circle (25+RadiusCirc,350),RadiusCirc,colorWhite
Circle (25+RadiusCirc+RadiusCirc+RadiusCirc+25,350),RadiusCirc,colorWhite

Locate 2 ,30 : Print "Location of Antennas"
Locate 2,105 : Print "Intensity in Direction"
Locate 20,2  : Print "^"
Locate 21,2  : Print "|"
locate 22,2  : print "N"
locate 42,42 : print "E->"

ExcentricityFactorX = 0.94

For n = 0 to  NumAnt   ' plot center and  antenna positions
    if n = 0 then
        x = 0 : y=0
    else
        x =  AntX(n) 
        y =  AntY(n) 
    end if
    x =  x*RadiusCirc/AntRadiusWave*ExcentricityFactorX + 25+RadiusCirc
    y = -y*RadiusCirc/AntRadiusWave + 350
    If n = 0 then 
        Line (X,Y+4)-(X,Y-4),colorRed
        Line (X+4,Y)-(X-4,Y),colorRed
      
    else
        Line (X,Y+2)-(X,Y-2),colorGreen
        Line (X+2,Y)-(X-2,Y),colorGreen
    End If
Next n



Locate 1,1         
Print "Number of antennas =";NumAnt, "Diameter in Wavelengths =";DiaWidthWave, "# of pairs =";NumAnt*(NumAnt-1)/2
 
Max = -1e20
Min = 1E20
 
 
DeltaRadians =  0.5/Rad2Deg  'Horizontal calib
DeltaRadians =  0.002*Deg2Rad  
DeltaRadians = RadiusField/RadiusCirc*Deg2Rad 
'DeltaRadians =  0.0002*Deg2Rad  
'DeltaRadians =  0.00002*Deg2Rad  
'DeltaRadians =  0.00002*Deg2Rad 
Locate 3,120
Print "Degrees/Pixel = ";DeltaRadians*Rad2Deg


Locate 44,90
Print "Normalized Power (P/n)"

y = 45 : x=90 : 
locate  y,x   : Color   colorWhite : print "1.0 - 0.3"
locate  y+1,x : Color   colorGray : print "0.3 - 0.1"
locate  y+2,x : Color   colorRed : print "0.1 - 0.03"
locate  y+3,x : Color   colorYellow : print "0.03 - 0.01"
locate  y+4,x : Color   colorGreen : print "0.01-0.003"
locate  y+5,x : Color   colorBlue : print "0.003-0.001"
locate  y+6,x : Color   colorViolet : print "0.001- "

Color colorWhite
 
 

color colorWhite
Locate 39,95 : Print "Magnified x3"
Locate 39,123 : Print "Magnified x3,"
Locate 40,123 : Print "Antenna beam width @ 2.3 GHz"
Locate  6,119 : Print "|<-------- 57 degrees -------->|"
Locate 38,89 : Print "|<------- 19 degrees -------->|"

' set up off axis average power 
TrigOffAxis = 0  
OffAxisPr = 0
OffAxisCnt = 0
OnAxisCnt = 0
'
LineCnt = 0    ' count lines printed
 
PhaseSum = 0
      HorOffsetPixels = 100 ' pixels
      LX = HorOffsetPixels
      LY = VertOffsetPixels+PhaseSum*10 
      PixelOffsetX = 0
      Locate   1,1
      PixelCnt = 0
For TargetXRadians = 0 to 1 step DeltaRadians ' outer angle loop
    PixelOffsetY = 0
    For TargetYRadians = 0 to 1 step DeltaRadians ' inner angle loop
       ' locate 54,54
     '   print "TargetXRadians=";TargetXRadians,"TargetYRadians";TargetYRadians;"           "  

        ' we want the normal of the wave plane that passes through the antenna !!!
        ' OK - we want three points of the plane to make the wave plane,
        '   we use the XZ plane angle (alpha) and the YZ plane angle (beta) 
        ' one point "P" is the antenna farm reference (0,0,0)
        ' nomenclature follows 
        ' http://www.math.washington.edu/~king/coursedir/m445w04/notes/vector/normals-planes.html#cross
        P1 = 0 : P2 = 0 : P3 = 0   ' using x,y,z
        '     - we have to cross product it cause Gauss/Jordan with the 0,0,0 point fails
        ' another point "X" is on the XZ plane, using sin (TargetXRadians)
        QXZ1x = cos(TargetXRadians) : QXZ2y = 0 : QXZ3z = sin(TargetXRadians)
        RYZ1x = 0 : RYZ2y = cos(TargetYRadians) : RYZ3z = sin(TargetYRadians)
        ' OK, now we have the three points in the plane wave front.
        '    now make two vectors 
        '  PQ = Q - P , but since P = 0, no change
        '  PR = R - P , but since P = 0, no change
        '
        ' OK, now the cross product to define the normal to the plane , the direction of the subject
        '     this is from Anton's "Elementary Linear Algebra, page 104
        ' U x V = (u2v3 - u3v2, u3v1 - u1v3,  u1v2 - u2v1)
        ' or changing names
        ' Q x R = (q2r3 - q3r2, q3r1 - q1r3,  q1r2 - q2r1)
        QxR1 = QXZ2y * RYZ3z - QXZ3z * RYZ2y
        QxR2 = QXZ3z * RYZ1x - QXZ1x * RYZ3z
        QxR3 = QXZ1x * RYZ2y - QXZ2y * RYZ1x
         
        'locate 55,55
   '     print "QxR1=";QxR1,"QxR2=";QxR2, "QxR3=";QxR3;"           " 
        ' OK, we have the normal to the plane   :-))
        '     the above 20  lines can obviously be optimised is several ways
        '          - but clarity seems more important
        ' now the final goal, to get the equation of the wave front, passing through 0,0,0
        '
        '     from page 113 on Anton (above), and point = 0,0,0 we get plane equation of
        '    QxR1*x + QxR2*y +QxR3*z + 0 = 0
        AbsoluteN = sqr(QxR1*QxR1 + QxR2*QxR2 + QxR3*QxR3)  ' not necessarily one, why???
       ' locate 56,56
    '    print "AbsoluteN =";AbsoluteN;"           " 
  '      input "input 'e' to stop";inkey$  ' this input requires a CR   :-((
  ' WHILE INKEY$ = "": sleep (10):WEND
   '     a$ = LCASE(inkey$) 
    '    print "I see ";a$
    ' ("e" <> LCASE(chr(getkey))  )   ' works
   '     if a$ = "e" then
  '        print    a$
  '        stop
  '      End if
        ' probably a reasonable optimization to do phasing to all antennas here,
        '      same repetitions later 
        LogWidthInt = Int(Log(DiaWidthWave))
    '    Print "LogWidthInt = ";LogWidthInt
          Locate   1,1
        For x = 1 to NumAnt   ' loop to set phase of every antenna
            ' Then with the wave front equation passing through 0,0,0,  
            '   from http://softsurfer.com/Archive/algorithm_0104/algorithm_0104.htm#Distance%20Point%20to%20Plane
            ' we get d(P0,pi) = n*(P0-V0)/|n| = (ax0 + by0 + cz0 + d)/srq(a^2 b^2+c^2)
            '   and if {n| = 1,
            '       we get d(P0,pi) = n*(P0-V0)/|n| = (ax0 + by0 + cz0 + d)
            ' for a plug in example - 
            ' from http://www.math.umn.edu/~nykamp/m2374/readings/planedistex/
            '  we get "P = (4, - 4, 3) to the plane 2x - 2y + 5z + 8 = 0.
            '           From the equation for the plane we substitute 
            '              A = 2, B = - 2, C = 5, D = 8. 
            '           From the point P, we substitute x1 = 4, y1 = - 4, and z1 = 3. 
            '    d = (|2*4 + (-2)*(-4) + 5*3 + 8|)/sqr(2^2 + (-2)^2 + 5^2) = roughly 6.8
        'QxR1 = QXZ2y * RYZ3z - QXZ3z * RYZ2y
        'QxR2 = QXZ3z * RYZ1x - QXZ1x * RYZ3z
        'QxR3 = QXZ1x * RYZ2y - QXZ2y * RYZ1x
        'AbsoluteN
            LogWidthInt = Int(Log(DiaWidthWave))
            DeltaWaveLengths = Abs((QxR1*AntX(x) + QxR2*AntY(x) + QxR3*0) / AbsoluteN)
            ' now we need the fractional part of the wavelength for phase 
            While  DeltaWaveLengths > LogWidthInt
                DeltaWaveLengths = DeltaWaveLengths - LogWidthInt
            WEnd
            While  DeltaWaveLengths >= 1
                DeltaWaveLengths = DeltaWaveLengths - 1
            WEnd
            AntPhase(x) = Cos(DeltaWaveLengths*TwoPi)
            'If (PixelOffsetX = 0) and (PixelOffsetY = 0) then 
          '      Print "AntPhase(x)=";AntPhase(x);"           "
           ' End If            
        Next x  ' end on set phase
       '  WHILE INKEY$ = "": sleep (10):WEND
      '  Print "AntPhase(NumAnt)=";AntPhase(NumAnt),"x = ";NumAnt;"           " 
        PairCnt = 0
        CorrelatedPairSum = 0 
        For BaseAnt = 1 to NumAnt    '  outer pair loop
           BasePhase = AntPhase(BaseAnt) 
           For PairAnt = BaseAnt+1 to NumAnt   ' inner pair loop
               CorrelatedPair =  BasePhase * AntPhase(PairAnt)     
               CorrelatedPairSum = CorrelatedPairSum + CorrelatedPair  ' integrating the voltage (power per antenna = 1 )
                ' Print "WkCos=";WkCos,"WkRad=";WkRad, "PhaseSum=";PhaseSum
               PairCnt = PairCnt+1
           Next PairAnt
    Next BaseAnt
    If CorrelatedPairSum < 0 then CorrelatedPairSum = 0
        
    PhaseSq = CorrelatedPairSum  / PairCnt   ' normalize to 1
    ' Print "PhaseSq =";PhaseSq;"   ";"CorrelatedPairSum";CorrelatedPairSum;"   ";"PairCnt";PairCnt;"      "
    If PhaseSq < 0.01 then TrigOffAxis = 1 
    If TrigOffAxis <> 0 then
        OffAxisPr = PhaseSq + OffAxisPr
        OffAxisCnt = OffAxisCnt + 1
    else
        OnAxisCnt = OnAxisCnt+1
    End if ' TrigOffAxis
    
    If PhaseSq > 0.3 then
        PixelColor = colorWhite
        PowerArrayByColor(1) = PowerArrayByColor(1) + 1
    ElseIf PhaseSq > 0.1 then
         PixelColor = colorGray
        PowerArrayByColor(2) = PowerArrayByColor(2) + 1
    ElseIf PhaseSq > 0.03 then
         PixelColor = colorRed
        PowerArrayByColor(3) = PowerArrayByColor(3) + 1
    ElseIf PhaseSq > 0.01 then
         PixelColor = colorYellow
        PowerArrayByColor(4) = PowerArrayByColor(4) + 1
    ElseIf PhaseSq > 0.003 then
         PixelColor = colorGreen
        PowerArrayByColor(5) = PowerArrayByColor(5) + 1
    ElseIf PhaseSq > 0.001 then
         PixelColor = colorBlue
        PowerArrayByColor(6) = PowerArrayByColor(6) + 1
    Else
        PixelColor = colorViolet
        PowerArrayByColor(7) = PowerArrayByColor(7) + 1
    End IF
   ' Print "PhaseSq=";PhaseSq;"              "
    PixelCenterX = 25+RadiusCirc+RadiusCirc+RadiusCirc+25 + PixelOffsetX
    PixelCenterY = 350 - PixelOffsetY
    Line     (PixelCenterX,PixelCenterY)-(PixelCenterX ,PixelCenterY ),PixelColor
    PixelCnt = PixelCnt + 1 
    
    ' now do a gain of 3
     if (PixelOffsetX < 80) and (PixelOffsetY < 80) then 
        PixelCenterX = 25+RadiusCirc+RadiusCirc+RadiusCirc+25 - PixelOffsetX*3
        PixelCenterY = 350 + PixelOffsetY*3
        Line     (PixelCenterX,PixelCenterY)-(PixelCenterX-3 ,PixelCenterY+3 ),PixelColor,BF
    End If

    ' now do antenna masking with gain 3
    if (PixelOffsetX < 80) and (PixelOffsetY < 80) then 
        ' estimated attenuation at 2.3 GHZ is 25R^2 in degrees with floor of 300 (about 25 db
        Atten = 10 * ( (TargetXRadians*Rad2Deg)^2 + (TargetYRadians*Rad2Deg)^2 )
        If Atten > 100 then 
            Atten = 100
        End If
        MaskPhaseSq = PhaseSq/Atten
        If MaskPhaseSq > 0.3 then
            PixelColor = colorWhite
            MaskPowerArrayByColor(1) = MaskPowerArrayByColor(1) + 1
        ElseIf MaskPhaseSq > 0.1 then
            PixelColor = colorGray
            MaskPowerArrayByColor(2) = MaskPowerArrayByColor(2) + 1
        ElseIf MaskPhaseSq > 0.03 then
            PixelColor = colorRed
            MaskPowerArrayByColor(3) = MaskPowerArrayByColor(3) + 1
        ElseIf MaskPhaseSq > 0.01 then
            PixelColor = colorYellow
            MaskPowerArrayByColor(4) = MaskPowerArrayByColor(4) + 1
        ElseIf MaskPhaseSq > 0.003 then
            PixelColor = colorGreen
            MaskPowerArrayByColor(5) = MaskPowerArrayByColor(5) + 1
        ElseIf MaskPhaseSq > 0.001 then
            PixelColor = colorBlue
            MaskPowerArrayByColor(6) = MaskPowerArrayByColor(6) + 1
        Else
            PixelColor = colorViolet
            MaskPowerArrayByColor(7) = MaskPowerArrayByColor(7) + 1
        End IF
        PixelCenterX = 25+RadiusCirc+RadiusCirc+RadiusCirc+25 + PixelOffsetX*3
        PixelCenterY = 350 + PixelOffsetY*3 + 1
        Line (PixelCenterX,PixelCenterY)-(PixelCenterX+3 ,PixelCenterY+3 ),PixelColor , BF 
    End If
  
        PixelOffsetY = PixelOffsetY + 1
    Next TargetYRadians   ' inner angle loop
    PixelOffsetX = PixelOffsetX + 1
Next TargetXRadians   ' outer angle loop


y = 45 : x=105 : 
locate  y,x   : Color   colorWhite :  print using "#.###^^^^";PowerArrayByColor(1)*100/PixelCnt
locate  y+1,x : Color   colorGray :   print  using "#.###^^^^";PowerArrayByColor(2)*100/PixelCnt
locate  y+2,x : Color   colorRed :    print  using "#.###^^^^";PowerArrayByColor(3)*100/PixelCnt
locate  y+3,x : Color   colorYellow : print  using "#.###^^^^";PowerArrayByColor(4)*100/PixelCnt
locate  y+4,x : Color   colorGreen :  print  using "#.###^^^^";PowerArrayByColor(5)*100/PixelCnt
locate  y+5,x : Color   colorBlue :   print  using "#.###^^^^";PowerArrayByColor(6)*100/PixelCnt
locate  y+6,x : Color   colorViolet : print  using "#.###^^^^";PowerArrayByColor(7)*100/PixelCnt        

Color colorWhite

' 
 
If OffAxisCnt > 0 then
    locate  y+7,90 
    Print  using "Off Axis average power #.###^^^^  for ###### values"; OffAxisPr/OffAxisCnt,OffAxisCnt 
 End if       
  WHILE INKEY$ = "": sleep (10):WEND
    STOP

WHILE INKEY$ = "": sleep (10):WEND
    STOP