pro thz_map,ptg=ptg,amp,pha,map,beam

  lambda = 3.333  ; Wavelength [mm]
  Z = 272.   ; Distance to source plane [mm]
  ; Positions of detectors [mm] (including baseline corrections)
  xo = [27.62,11.9,-15.98,11.9] + [1.4,0.7,-6.1,0.7]
  yo = [0,-37.12,0,20.93] + [1.4,-1.5,4.5,15.9]
  ; Nominal positions of source
  x = [0,-10,10,10,-10]         ; X position of source for 5 different positions
  y = [0,-10,-10,10,10]         ; Y position of source for 5 different positions

  if (not keyword_set(ptg)) then ptg = 0
  ; Correct phase of each detector for near-field
  corpha = pha
  ; Loop over source positions
  ;for j = 0, 4 do begin
     ; Loop over the four detectors
     for i = 0, 3 do begin
        ; Angle of detector i from source
;        theta = atan(sqrt((xo[i]-x[ptg])^2+(yo[i]-y[ptg])^2),Z)
        theta = atan(sqrt((xo[i])^2+(yo[i])^2),Z)  ; It seems to be wrong to correct for the source posn shifts
        delphi = 2*!dpi*(Z/lambda)*(1/cos(theta) - 1)  ; Phase shift at detector i
        case ptg of
           0: corpha[i,1550:2350] = pha[i,1550:2350]-delphi
           1: corpha[i,3200:4000] = pha[i,3200:4000]-delphi
           2: corpha[i,4080:4750] = pha[i,4080:4750]-delphi
           3: corpha[i,4870:5550] = pha[i,4870:5550]-delphi
           4: corpha[i,5670:6249] = pha[i,5670:6249]-delphi
        endcase
     endfor
  ;endfor
  
  ; Declare storage for u,v and data
  nbl = 6   ; Number of baselines for 4 detectors
  nt = n_elements(amp[0,*])
  uvdat = complexarr(nbl,n_elements(amp[0,*]))
  u = (v = fltarr(nbl)) 
  bp = fltarr(nbl)   ; To hold beam phase near-field corrections
  ; Calculate baseline amplitudes and phases, using corrected phase
  k = -1  ; Baseline index
  for i = 0, 2 do begin
     for j = i+1, 3 do begin
        k = k + 1
        u[k] = (xo[j]-xo[i])/lambda
        v[k] = (yo[j]-yo[i])/lambda
        a = sqrt(amp[i,*]*amp[j,*])
        p = corpha[j,*]-corpha[i,*]
        uvdat[k,*] = complex(a*cos(p),a*sin(p))
     endfor
  endfor
  
  ; Generate image using discrete FT rather than FFT
  N = 100
  xshift = 0
  yshift = 0
  xi=(findgen(N+1)-N/2+xshift)#(fltarr(N+1)+1)/Z
  eta=transpose((findgen(N+1)-N/2+yshift)#(fltarr(N+1)+1))/Z
  xi2 = (findgen(2*N+1)-N+xshift)#(fltarr(2*N+1)+1)/Z
  eta2 = transpose((findgen(2*N+1)-N+yshift)#(fltarr(2*N+1)+1))/Z
  map = fltarr(N+1,N+1)
  beam = fltarr(2*N+1,2*N+1)
  case ptg of
    0: begin
         for it = 1550, 2350 do begin
            for k = 0, nbl-1 do begin
               map = map + float(uvdat[k,it])*cos(2*!pi*(u[k]*xi+v[k]*eta))-imaginary(uvdat[k,it])*sin(2*!pi*(u[k]*xi+v[k]*eta))
               beam = beam + 1*cos(2*!pi*(u[k]*xi2+v[k]*eta2))
            endfor
         endfor
       end
    1: begin
         for it = 3200, 4000 do begin
            for k = 0, nbl-1 do begin
               map = map + float(uvdat[k,it])*cos(2*!pi*(u[k]*xi+v[k]*eta))-imaginary(uvdat[k,it])*sin(2*!pi*(u[k]*xi+v[k]*eta))
               beam = beam + 1.*cos(2*!pi*(u[k]*xi2+v[k]*eta2))
            endfor
         endfor
       end
    2: begin
         for it = 4080, 4750 do begin
            for k = 0, nbl-1 do begin
               map = map + float(uvdat[k,it])*cos(2*!pi*(u[k]*xi+v[k]*eta))-imaginary(uvdat[k,it])*sin(2*!pi*(u[k]*xi+v[k]*eta))
               beam = beam + 1.*cos(2*!pi*(u[k]*xi2+v[k]*eta2))
            endfor
         endfor
       end
    3: begin
         for it = 4870, 5550 do begin
            for k = 0, nbl-1 do begin
               map = map + float(uvdat[k,it])*cos(2*!pi*(u[k]*xi+v[k]*eta))-imaginary(uvdat[k,it])*sin(2*!pi*(u[k]*xi+v[k]*eta))
               beam = beam + 1.*cos(2*!pi*(u[k]*xi2+v[k]*eta2))
            endfor
         endfor
       end
    4: begin
         for it = 5670, 6249 do begin
            for k = 0, nbl-1 do begin
               map = map + float(uvdat[k,it])*cos(2*!pi*(u[k]*xi+v[k]*eta))-imaginary(uvdat[k,it])*sin(2*!pi*(u[k]*xi+v[k]*eta))
               beam = beam + 1.*cos(2*!pi*(u[k]*xi2+v[k]*eta2))
            endfor
         endfor
       end
  endcase       
  
return
end