;-------- wind speed in m/s ------------------------------------
function fwind, altitude, latitude, dat_wind, ln_wind

if altitude ge dat_wind[0,ln_wind-1] then begin
  vw = (dat_wind[1,ln_wind-1] - dat_wind[1,ln_wind-2]) / (dat_wind[0,ln_wind-1] - dat_wind[0,ln_wind-2]) * (altitude - dat_wind[0,ln_wind-1]) + dat_wind[1,ln_wind-1]
  vw = vw * cos(!dpi * latitude / 180.) / cos(!dpi * 45. / 180.)
  return, vw
endif
if altitude le 0.0 then return, 0.0
i = 0
while dat_wind[0,i] lt altitude do i = i + 1
vw = (dat_wind[1,i] - dat_wind[1,i-1]) / (dat_wind[0,i] - dat_wind[0,i-1]) * (altitude - dat_wind[0,i-1]) + dat_wind[1,i-1]
vw = vw * cos(!dpi * latitude / 180.) / cos(!dpi * 45. / 180.)
return, vw
end
;----- end of fwind --------------------------------------------

;-------- Cassini position ------------------------------------
function fcpos, time, dat_cas, ln_cas, index

if time ge dat_cas[0,ln_cas-1] then begin
  alt = (dat_cas[index,ln_cas-1] - dat_cas[index,ln_cas-2]) / (dat_cas[0,ln_cas-1] - dat_cas[0,ln_cas-2]) * (time - dat_cas[0,ln_cas-1]) + dat_cas[index,ln_cas-1]
  return, alt
endif
;if time le 0.0 then return, 0.0
i = 0
while dat_cas[0,i] lt time do i = i + 1
alt = (dat_cas[index,i] - dat_cas[index,i-1]) / (dat_cas[0,i] - dat_cas[0,i-1]) * (time - dat_cas[0,i-1]) + dat_cas[index,i-1]
return, alt
end
;----- end of fcpos --------------------------------------------

;----- start of main program --------------------------------------------
print, 'Output Medium:                      '
print, '                                    '
print, '               output.ps         (0)'
print, '               output.eps        (1)'
print, '               Tektronix-Screen  (2)'
print, '               X-Screen          (3)'
print, '               Windows           (4)'
print, '               terminate program (5)'
read, medium

case medium of
  0: begin
       set_plot, 'PS'
       device, filename='output.ps', $
         bits=8, /color, $
         /landscape,  $
         xoffset=2.5, yoffset=26., $
         xsize=24., ysize=16.

;         /portrait, $
;         xoffset=2., yoffset=3., $
;         xsize=16., ysize=24.
     end
  1: begin
       set_plot, 'PS'
       device, filename='output.eps', $
         bits=8, /color, $
         /encapsulated, $
         xsize=24., ysize=16.
     end
  2: begin
       set_plot, 'TEK'
     end
  3: begin
       set_plot, 'X'
;       window, 0, xsize=650.*210/297, ysize=650.
;       !p.color=255
;       !p.background=0
     end
  4: begin
       set_plot, 'WIN'
;       window, 0, xsize=600., ysize=400.
;       !p.color=255
;       !p.background=0
     end
  else: begin
          print, 'terminating program'
          stop
        end
endcase

xyouts, 0.5, 0.5, /normal, '!5'
dummy = ' '

RED   = [0, 1, 1, 0, 0, 1, 1]
GREEN = [0, 1, 0, 1, 0, 1, 0.5]
BLUE  = [0, 1, 0, 0, 1, 0, 0]
tvlct, 255 * RED, 255 * GREEN, 255 * BLUE
black=0 & white=1 & red=2 & green=3 & blue=4 & yellow=5 & orange=6

;-------- landing site co-ordinates ---------
tour = 1
;read, tour, prompt = 'Nominal (1) or Backup Tour (2)? '
if tour eq 1 then begin
  lon_ls0 = 152.0  ; west
  lat_ls0 = 18.75  ; north
  print, 'using landing coordinates of nominal tour.'
endif else begin
  lon_ls0 = 146.26 ; west
  lat_ls0 = 24.43  ; north
  print, 'using landing coordinates of back-up tour.'
endelse
print, 'Huygens landing site:'
print, 'latitude : ' + string(lat_ls0) + ' deg'
print, 'longitude: ' + string(lon_ls0) + ' deg'
;-------- done: landing site co-ordinates ---------

;-------- inputs ----------
read, t_eval, prompt = 'Enter evaluation time (min after t0): '

; Cassini's pointing on Titan's surface
;read, caspoint_lat, prompt = 'Enter latitude for Cassini HGA pointing (deg): '
;read, caspoint_lon, prompt = 'Enter longitude for Cassini HGA pointing (deg): '
caspoint_lat = lat_ls0
caspoint_lon = lon_ls0
;caspoint_lon = 140.
print, 'Cassini pointing to:'
print, 'latitude : ' + string(caspoint_lat) + ' deg'
print, 'longitude: ' + string(caspoint_lon) + ' deg'
;-------- done: inputs ----------

;------- read wind model -----------------------------
filename = 'flasar.dat'
;read, filename, prompt = 'Enter filename for wind model: '
;openr, 1, filename
;ln_wind = 0
;while not eof(1) do begin
;  readf, 1, dummy
;  ln_wind = ln_wind + 1
;endwhile
;close, 1
ln_wind = 20
dat_wind = dblarr(2,ln_wind)
openr, 1, filename
readf, 1, dat_wind
close, 1
;------ done: read wind model -----------------------------

;------- read HUYGENS trajectory data -----------------------------
filename = 'nom_des.dat'
;read, filename, prompt = 'Enter filename for HUYGENS trajectory data: '
;openr, 1, filename
;ln_huy = 0
;while not eof(1) do begin
;  readf, 1, dummy
;  ln_huy = ln_huy + 1
;endwhile
;close, 1
ln_huy = 482
dat_huy = dblarr(13,ln_huy)
openr, 1, filename
readf, 1, dat_huy
close, 1
; index for dat_huy:
; 0: time (t-t0 min), 1: altitude (km), 2: vert. vel. (m/s), 3: acc (m/s2),
; 4: range (km), 5: r-r (m/s), 6: alfa, 7: prl, 8: xlat,
; 9: paa (PAA = paa + alfa), 10: sza, 11: dwc, 12: phi
;------- done: read HUYGENS trajectory data -----------------------------

;------- read CASSINI trajectory data -----------------------------
filename = 'orbiter-flyby_t0.dat'
;read, filename, prompt = 'Enter filename for CASSINI trajectory data: '
;openr, 1, filename
;ln_cas = 0
;while not eof(1) do begin
;  readf, 1, dummy
;  ln_cas = ln_cas + 1
;endwhile
;close, 1
ln_cas = 50
dat_cas = dblarr(13,ln_cas)
openr, 1, filename
readf, 1, dat_cas
close, 1
; index for dat_cas:
; 0: time (t-t0 min), 1: range (km), 2: altitude (km), 3: ssc_lat,
; 4: ssc_lon, 5: ss_lat, 6: ss_lon, 7: v_rad, 8: v_tan, 9: phase,
; 10: AD (rad), 11: R.A., 12: dec
;------- done: read CASSINI trajectory data -----------------------------

;-------- constants ---------------------------------
RT = 2575.0             ; Titan radius (km)
sigma_NS = 50.0 / RT   ; 1 sigma delivery error ellipse, NS axis (rad)
sigma_EW = 160.0 / RT / cos(lat_ls0 * !dpi/180.) ; 1 sigma delivery error ellipse, EW axis (rad)
no_bases = 1000         ; number of bases for Monte Carlo
;-------- end: constants ---------------------------------

;-------- initial error ellipse ---------------------------------
lat_ls = randomn(seed, no_bases) * sigma_NS * 180.0/!dpi + lat_ls0 ; (deg)
lon_ls = randomn(seed, no_bases) * sigma_EW * 180.0/!dpi + lon_ls0 ; (deg)
;-------- done: initial error ellipse ---------------------------------

;-------- evolution of error ellipse ---------------------------------
if t_eval ge dat_huy[0,ln_huy-1] then begin
  t_final = dat_huy[0,ln_huy-1]  ; moment for evaluation of error ellipse
  no_ttags = ln_huy        ; number of time tags for wind integration
endif else begin
  i = 0
  while dat_huy[0,i] lt t_eval do i = i + 1
  t_final = dat_huy[0,i]
  no_ttags = i + 1
endelse

v_wind = dblarr(no_bases, no_ttags)
;lat    = dblarr(no_bases, no_ttags)  ; no change in lat with time !!
lon    = dblarr(no_bases, no_ttags)

v_randomn = randomn(seed, no_bases)  ; random number for wind speed generation

;-------------- manipulate v_wind here !!----------------
scfl = 1.5  ; positive <=> prograde, negative <=> retrograde
sigma_wind = 0.5
retro = 0 ; 1: 6% retrograde winds, 0: 100% prograde

if scfl ge 0. then begin
  if retro then dummy = '!C(94% prograde, 6% retrograde)' $
  else dummy = '!C(100% prograde)'
  form = '(d3.1)'
endif else begin
  if retro then dummy = '!C(94% retrograde, 6% prograde)' $
  else dummy = '!C(100% retrograde)'
  form = '(d4.1)'
endelse

wms = 'Wind model: (' + string(scfl, format = form) + ' !9+!X ' + $
  string(sigma_wind, format = '(d4.2)') + ') * Flasar-Wind' + dummy
;-------------- done: manipulate v_wind here !!----------------

for i = 0, no_bases - 1 do begin
  for j = 0, no_ttags - 1 do begin

    v_wind[i,j] = fwind(dat_huy[1,j], lat_ls[i], dat_wind, ln_wind) * (scfl + sigma_wind * v_randomn[i])

    if retro and ((i + 1) gt (0.94 * no_bases)) then $   ; 6% probability for retrograde winds
      v_wind[i,j] = -v_wind[i,j]

; v_wind positive <=> eastward, but delta_lon positive <=> westward !!
; v_wind in m/sec !!

    if j ge 1 then $
      lon[i,j] = lon[i,j-1] - 27./!dpi/5. * (v_wind[i,j-1] / (RT + dat_huy[1,j-1]) + v_wind[i,j] / (RT + dat_huy[1,j])) * (dat_huy[0,j] - dat_huy[0,j-1]) $
    else lon[i,j] = lon_ls[i]

  endfor
endfor

;-------- plot wind model -----------------------

!p.multi = [0, 1, 2, 0, 0]

plot, abs(scfl) * dat_wind[1,*] * cos(!dpi * lat_ls0 / 180.) / 0.707106781187, dat_wind[0,*], $
;plot, dat_wind[1,*] * cos(!dpi * lat_ls0 / 180.) / cos(!dpi * 45. / 180.), dat_wind[0,*], $
  title = wms, $
  xtitle = 'MEAN ABSOLUTE WIND VELOCITY [m/sec]', $
  ytitle = 'ALTITUDE [km]', $
  xstyle = 1, $
;  yrange = [0., 170.], $
  ystyle = 1, $
  color = black

xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/20., $
        !y.crange[0]+16.*(!y.crange[1]-!y.crange[0])/20., $
  'Latitude: ' + string(lat_ls0, format = '(d5.2)') + '!9%!XN', $
  color = black

bins = 10
xx = indgen(ceil((max(v_wind[*,0]) - min(v_wind[*,0])) / bins)) * bins
xx = xx + floor(min(v_wind[*,0]))
yy = 100. * histogram(v_wind[*,0], binsize = bins) / no_bases

;plot, xx, 100. * histogram(v_wind[*,0], binsize = bins) / no_bases, $
plot, xx, yy, $
  psym = 10, $
  title = 'Wind Velocity Density Function at ' + string(dat_huy[1,0], format = '(d5.1)') + ' km Altitude', $
  xtitle = 'WIND VELOCITY [m/sec]', $
  ytitle = 'PERCENTAGE (BINSIZE = ' + string(bins, format = '(i2)') + ' m/s)', $
  xstyle = 1, $
;  ystyle = 1, $
;  ystyle = 3, $
  color = black
;-------- done: plot wind model -----------------------
!p.multi = 0

xmin = 200 & xmax = 100 & ymin = -10 & ymax = 50

plot, lon[*,no_ttags-1], lat_ls[*], /nodata, color=black, $
  title = "Cassini-HGA Beam (2040 MHz) on Titan's Surface at !18t-t!D0!N!X = " + string(t_eval, format='(i3)') + " min", $
;  title = "Cassini-HGA Beam (2098 MHz) on Titan's Surface at !18t-t!D0!N!X = " + string(t_eval, format='(i3)') + " min", $
  xstyle = 1, ystyle = 1, $
  ytitle = 'LATITUDE [!9%!XN]', $
  xtitle = 'LONGITUDE [!9%!XW]', $
  xrange = [xmin, xmax], $
  yrange = [ymin, ymax]
;oplot, lon[*,0], lat_ls[*], psym=3, color=blue
oplot, lon[*,no_ttags-1], lat_ls[*], psym=3, color=green

xyouts, xmin+(xmax-xmin)/20., ymin+18.*(ymax-ymin)/20., $
  'Huygens at !18t!D0!N!X (Nominal: ' + string(lat_ls0, format='(d5.2)') + $
  '!9%!XN, ' + string(lon_ls0, format='(d6.2)') + '!9%!XW)', $
  color = blue
xyouts, xmin+(xmax-xmin)/20., ymin+16.*(ymax-ymin)/20., $
  'Huygens at !18t-t!D0!N!X = ' + string(t_eval, format='(i3)') + ' min', $
  color = green

;-------- done: evolution of error ellipse ---------------------------------

;-------- projection of antenna pattern on Titan's surface -----------------
caspos = dblarr(3)
caspoint_pos = dblarr(3)
c2h = dblarr(3)
absc2h = dblarr(3)   ; need 2 arrays for determination of r-r, and 1 for t0-value

delta_t = 2. ; minutes

for i = 0, 2 do begin       ; nec. for r-r and t0-value

  case i of
    0: begin
       tt = 0
       end
    1: begin
       tt = t_eval + delta_t
       end
    2: begin
       tt = t_eval
       end
  endcase

  ; Cassini position at t_eval <=> i = 2
  ;casalt = fcpos(tt, dat_cas, ln_cas, 2)
  casalt = fcpos(tt, dat_cas, ln_cas, 1)     ; range, not altitude !!
  caslat = fcpos(tt, dat_cas, ln_cas, 3)
  caslon = fcpos(tt, dat_cas, ln_cas, 4)

  sintheta = sin(!dpi/2. - caslat * !dpi/180.)     ; theta = 0 deg: north pole (lat = 90 deg)
  costheta = cos(!dpi/2. - caslat * !dpi/180.)     ; theta = 180 deg: south pole (lat = -90 deg)
  sinphi   = 0.                                    ; phi = 0 deg: longitude of Cassini
  cosphi   = 1.

  ;Cassini position, cartesian
  caspos[0] = casalt * sintheta * cosphi
  caspos[1] = casalt * sintheta * sinphi
  caspos[2] = casalt * costheta

  sintheta = sin(!dpi/2. - caspoint_lat * !dpi/180.)
  costheta = cos(!dpi/2. - caspoint_lat * !dpi/180.)
  sinphi   = sin((caspoint_lon - caslon) * !dpi/180.)
  cosphi   = cos((caspoint_lon - caslon) * !dpi/180.)

  ;Cassini pointing on Titan's surface, cartesian
  caspoint_pos[0] = RT * sintheta * cosphi
  caspoint_pos[1] = RT * sintheta * sinphi
  caspoint_pos[2] = RT * costheta

  ;Cassini pointing direction, cartesian vector Cassini -> Huygens
  c2h = caspoint_pos - caspos
  absc2h[i] = sqrt(c2h[0]*c2h[0] + c2h[1]*c2h[1] + c2h[2]*c2h[2])
endfor

c2h_dot = (absc2h[1] - absc2h[2]) / delta_t / 60.
; r-r in km/sec, note: Huygens is not moving!!

;--------- assign power value to each Huygens --------------------

huy_pos = dblarr(3, no_bases)
huy_offang = dblarr(no_bases)
huy_powdb = dblarr(no_bases)

for i = 0, no_bases - 1 do begin
  sintheta = sin(!dpi/2. - lat_ls[i] * !dpi/180.)
  costheta = cos(!dpi/2. - lat_ls[i] * !dpi/180.)
  sinphi   = sin((lon[i, no_ttags-1] - caslon) * !dpi/180.)
  cosphi   = cos((lon[i, no_ttags-1] - caslon) * !dpi/180.)

  ;Huygens position on Titan's surface, cartesian
  huy_pos[0,i] = RT * sintheta * cosphi
  huy_pos[1,i] = RT * sintheta * sinphi
  huy_pos[2,i] = RT * costheta

  cosx = total(c2h[*] / RT * (huy_pos[*,i] - caspos[*]) / RT) / $
     sqrt(total(c2h[*] / RT * c2h[*] / RT) * $
     total((huy_pos[*,i] - caspos[*]) / RT * $
           (huy_pos[*,i] - caspos[*]) / RT))
  if (cosx gt 1.) or (cosx lt (-1.)) then begin
    print, 'error while assigning power value for each Huygens:'
    print, 'i = ' + string(i)
    print, 'cos(x) = '+ string(cosx)
    print, ' '
    huy_offang[i] = 0.0
  endif $
  else huy_offang[i] = acos(cosx) * 180./!dpi

endfor

;;huy_powdb[i] = 10 * log 2 * x^2 / x_3db^2    ; 3dB: 1.2 deg
;huy_powdb = 3.0103 * huy_offang * huy_offang / 1.44
; s.u.

;--------- done: assign power value to each Huygens --------------------

;ak = dblarr(2)
;ak[0] = 1.2  ; half 3-dB-beam (deg) for chain A
;ak[1] = 1.67  ; half 6-dB-beam (deg) for chain A
;ak[0] = 1.14  ; half 3-dB-beam (deg) for chain B
;ak[1] = 1.56  ; half 6-dB-beam (deg) for chain B
; a-keule (ak) defined via a-power (ap),
; ap follows gauss-function, see for-loop

;no_keulen = 6
no_keulen = 8
ak = dblarr(no_keulen)  ; for 3, 6, 9, 12, 15, 18 dB beams, deg
;ap = dblarr(no_keulen)  ; for 3, 6, 9, 12, 15, 18 dB beams, inverse power

;no_kv = 60
no_kv = 30

a = dblarr(3, no_kv)
; a-vector for construction of beam, perp. to c2h
; beam vectors are c2h+a
; no_kv vectors for psi = 0, 360, increment 360/no_kv
; initializiation doesn't work, if beam is parralel to y-axis (eq.-plane)!!

; calculation of lambda = extension or reduction factor for c2h+a
; condition: abs( caspos + lambda * (c2h + a) ) = R
l = dblarr(2)
kp = dblarr(3, no_kv)
lat_kp = dblarr(no_kv)
lon_kp = dblarr(no_kv)
lambda = dblarr(no_kv)

for keule = 0, no_keulen-1 do begin

;  ap[keule] = 2.^(keule + 1)                     ; inverse power, s.a.

  ak[keule] = sqrt(keule+1) * 1.2   ; 3-dB-beam for chain A: 1.20 deg
;  ak[keule] = sqrt(keule+1) * 1.14   ; 3-dB-beam for chain B: 1.14 deg

  tanak = tan(ak[keule] * !dpi/180.)

  absa = tanak * absc2h[2]
  a[2,0] = 0.
  a[1,0] = absa / sqrt(1.+c2h[1]/c2h[0]*c2h[1]/c2h[0])
  a[0,0] = a[1,0] * (-c2h[1])/c2h[0]

  ; construction of array set a (perp. to c2h)
  coef_a = ((-c2h[1])*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0])/(c2h[0]*(a[0,0]*(-c2h[1])/c2h[2]+a[1,0]))+c2h[2]/c2h[0])* $
    ((-c2h[1])*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0])/(c2h[0]*(a[0,0]*(-c2h[1])/c2h[2]+a[1,0]))+c2h[2]/c2h[0])+ $
    (a[0,0]*(-c2h[2])/c2h[0]+a[2,0])/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])* $
    (a[0,0]*(-c2h[2])/c2h[0]+a[2,0])/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])+1.

  for i = 1, no_kv/2 do begin

    psi = 360. / double(no_kv) * i
    cospsi = cos(psi * !dpi/180.)

; coef_a calculated outside for-loop !!

    coef_b = -(2.*(-c2h[1])*cospsi*absa*absa/(c2h[0]*(a[0,0]*(-c2h[1])/c2h[2]+a[1,0]))* $
      ((-c2h[1])*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0])/(c2h[0]*(a[0,0]*(-c2h[1])/c2h[2]+a[1,0]))+c2h[2]/c2h[0])+ $
      2.*cospsi*absa*absa*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0])/ $
      (a[0,0]*(-c2h[1])/c2h[2]+a[1,0])/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0]))

    coef_c = cospsi*absa*absa/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])* $
      cospsi*absa*absa/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])* $
      (1.+c2h[1]*c2h[1]/c2h[0]/c2h[0])-absa*absa

    sqrt_arg = coef_b*coef_b - 4.*coef_a*coef_c
    if sqrt_arg lt 0.0 then begin
      print, 'Construction of a:'
      print, 'sqrt_arg < 0 in a-calculation => a imaginary'
      print, 'i = ' + string(i) + ' (if i == no_kv/2 then expected)'
      print, 'no_kv = ' + string(no_kv)
      print, 'sqrt_arg = ' + string(sqrt_arg)
      if i eq no_kv / 2 then begin
        print, 'setting a[*,i] to -a[*,0]...'
        print, ' '
        a[2,i] = 0.0
        a[1,i] = -a[1,0]
        a[0,i] = -a[0,0]
      endif else begin
        print, 'setting sqrt_arg 0...'
        print, ' '
;        sqrt_arg = 0.
        a[2,i] = (-coef_b) / 2. / coef_a
        a[1,i] = (cospsi*absa*absa-a[2,i]*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0]))/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])
        a[0,i] = ((-c2h[1])*a[1,i] - c2h[2]*a[2,i]) / c2h[0]
      endelse
    endif else begin
      a[2,i] = (-coef_b + sqrt(sqrt_arg)) / 2. / coef_a
      a[1,i] = (cospsi*absa*absa-a[2,i]*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0]))/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])
      a[0,i] = ((-c2h[1])*a[1,i] - c2h[2]*a[2,i]) / c2h[0]
    endelse
    if i lt no_kv/2 then begin
      a[2,no_kv-i] = (-coef_b - sqrt(sqrt_arg)) / 2. / coef_a
      a[1,no_kv-i] = (cospsi*absa*absa-a[2,no_kv-i]*(a[0,0]*(-c2h[2])/c2h[0]+a[2,0]))/(a[0,0]*(-c2h[1])/c2h[2]+a[1,0])
      a[0,no_kv-i] = ((-c2h[1])*a[1,no_kv-i] - c2h[2]*a[2,no_kv-i]) / c2h[0]
    endif

  endfor

  for i = 0, no_kv-1 do begin

    if finite(a[2,i]) then begin
;      coef_a = (c2h[0] + a[0,i]) * (c2h[0] + a[0,i]) + $
;        (c2h[1] + a[1,i]) * (c2h[1] + a[1,i]) + (c2h[2] + a[2,i]) * (c2h[2] + a[2,i])
;      coef_b = 2. * (caspos[0] * (c2h[0] + a[0,i]) + caspos[1] * (c2h[1] + a[1,i]) + $
;       caspos[2]*(c2h[2]+a[2,i]))
;      coef_c = caspos[0]*caspos[0] + caspos[1]*caspos[1] + caspos[2]*caspos[2] - RT*RT

; coef_x too large for determination of lambda => division by RT*RT
      coef_a = ((c2h[0] + a[0,i]) * (c2h[0] + a[0,i]) + $
        (c2h[1] + a[1,i]) * (c2h[1] + a[1,i]) + (c2h[2] + a[2,i]) * (c2h[2] + a[2,i])) / RT/RT
      coef_b = (2. * (caspos[0] * (c2h[0] + a[0,i]) + caspos[1] * (c2h[1] + a[1,i]) + $
        caspos[2]*(c2h[2]+a[2,i]))) / RT/RT
      coef_c = (caspos[0]*caspos[0] + caspos[1]*caspos[1] + caspos[2]*caspos[2]) / RT/RT - 1.

      sqrt_arg = coef_b*coef_b - 4.*coef_a*coef_c

      if sqrt_arg ge 0. then begin
        l[0] = (-coef_b + sqrt(sqrt_arg)) / 2. / coef_a
        l[1] = (-coef_b - sqrt(sqrt_arg)) / 2. / coef_a
        if finite(l[0]) then lambda[i] = min(l) else begin
          lambda[i] = 1.e11
          print, 'Construction of lambda:'
          print, 'i = ' + string(i)
          print, 'error in lambda calculaton:'
          print, 'lambda infinite, because'
  	  print, 'coef_a = ' + string(coef_a)
          print, ' '
        endelse
      endif else begin
        lambda[i] = 1.e11
        print, 'Construction of lambda:'
        print, 'i = ' + string(i)
        print, 'lambda imaginary => beam array does not hit Titan'
        print, ' '
      endelse
    endif else begin
      lambda[i] = 1.e12
      print, 'Construction of lambda:'
      print, 'i = ' + string(i)
      print, 'a-array imaginary <= what does that mean???'
      print, ' '
    endelse

;    for j = 0, 2 do kp[j,i] = caspos[j] + lambda[i] * (c2h[j] + a[j,i])
    kp[*,i] = caspos[*] + lambda[i] * (c2h[*] + a[*,i])

    if abs(kp[2,i]/RT) le 1. then begin
      theta_kp = 180./!dpi * acos(kp[2,i]/RT)
      lat_kp[i] = 90.0 - theta_kp
    endif else begin
      print, 'Construction of latitude:'
      print, 'i = ' + string(i)
      print, 'kp[2,i] / RT = ' + string(kp[2,i] / RT)
      print, ' '
      theta_kp = 1.e10
      lat_kp[i] = 1.e10
    endelse

;  sqrt_arg = (RT^2 - (kp[2,i]^3)/RT) / (kp[0,i]^2 + kp[1,i]^2 * (2. + (kp[1,i]/kp[0,i])^2))

    if theta_kp le 180.0 then $
      lon_kp[i] = 180.0/!dpi * atan(kp[1,i], kp[0,i]) + caslon $
    else begin
      print, 'Construction of longitude:'
      print, 'i = ' + string(i)
      print, 'error in calculation of polar beam coordinates!'
      print, 'theta_kp = ' + string(theta_kp)
      print, ' '
      lon_kp[i] = 1.e10
      lat_kp[i] = 1.e10
    endelse
  endfor
;-------- done: projection of antenna pattern on Titan's surface -----------------

  oplot, lon_kp, lat_kp, max_value = 1000.0, $
    color = red, $
    linestyle = (keule lt 6) * keule + (keule ge 6) * (keule - 5)
  oplot, [lon_kp[no_kv-1], lon_kp[0]], [lat_kp[no_kv-1], lat_kp[0]], max_value = 1000.0, $
    color = red, $
    linestyle = (keule lt 6) * keule + (keule ge 6) * (keule - 5)

  if (lat_kp[no_kv/4] gt ymin) and (lat_kp[no_kv/4] lt ymax) and $
     (lon_kp[no_kv/4] gt xmax) and (lon_kp[no_kv/4] lt xmin) then $
    xyouts, lon_kp[no_kv/4], lat_kp[no_kv/4], $ 
; /4: oben, 3/4: unten, 0: links, /2: rechts
      string((keule + 1) * 3, format = '(i2)') + ' dB', $
      color = orange, charsize = 0.7, alignment = 0.5 $
;  else if (lat_kp[2*no_kv/3] gt ymin) and (lat_kp[2*no_kv/3] lt ymax) and $
;     (lon_kp[2*no_kv/3] gt xmax) and (lon_kp[2*no_kv/3] lt xmin) then $
;    xyouts, lon_kp[2*no_kv/3], lat_kp[2*no_kv/3], $ 
;      string((keule + 1) * 3, format = '(i2)') + ' dB', $
;      color = orange, charsize = 0.7, alignment = 0.5 $
  else if (lat_kp[no_kv/2+5] gt ymin) and (lat_kp[no_kv/2+5] lt ymax) and $
     (lon_kp[no_kv/2+5] gt xmax) and (lon_kp[no_kv/2+5] lt xmin) then $
    xyouts, lon_kp[no_kv/2+5], lat_kp[no_kv/2+5], $ 
      string((keule + 1) * 3, format = '(i2)') + ' dB', $
      color = orange, charsize = 0.7, alignment = 0.5

; always, instead of blue dots
;  if t_eval lt 0.9 then begin
    t = dindgen(401) * 2.0 * !DPI / 400.0
    x = 3. * sigma_EW * 180./!dpi * cos(t)
    y = 3. * sigma_NS * 180./!dpi * sin(t)
    oplot, lon_ls0 + x, lat_ls0 + y, $
      color = blue
;  endif

;-------- calculation of probabilities -----------------
;  count = no_bases
  tmp = where(huy_offang le ak[keule], count)
;
;  for i = 0, no_bases - 1 do begin  ; for each Huygens
;    sign_initial = 0
;    for j = 1, no_kv - 1 do begin
;      while (j lt no_kv - 1) and ((lon_kp[j-1] gt 1.e9) or (lon_kp[j] gt 1.e9) $
;        or (lat_kp[j-1] gt 1.e9) or (lat_kp[j] gt 1.e9)) $
;	do j = j + 1
;      x1 = lon_kp[j-1] - lon[i, no_ttags-1]
;      x2 = lon_kp[j] - lon[i, no_ttags-1]
;      y1 = lat_kp[j-1] - lat_ls[i]
;      y2 = lat_kp[j] - lat_ls[i]
;      if (x1 * y2 - y1 * x2) gt 0. then sign_curr = 1 else sign_curr = -1
;;      if j eq 1 then sign_initial = sign_curr $
;      if sign_initial eq 0 then sign_initial = sign_curr $
;      else if sign_curr ne sign_initial then begin
;        count = count - 1
;        j = no_kv - 1
;      endif
;    endfor
;  endfor

  prob = double(count) / double(no_bases)

  if (lat_kp[3*no_kv/4] gt ymin) and (lat_kp[3*no_kv/4] lt ymax) then $
    xyouts, lon_kp[3*no_kv/4], lat_kp[3*no_kv/4], $
; /4: oben, 3/4: unten, 0: links, /2: rechts
      string(prob * 100., format = '(d5.1)') + '%', $
      color = green, charsize = 0.7, alignment = 0.5

;-------- done: calculation of probabilities -----------------

endfor  ; keule

xyouts, xmin+(xmax-xmin)/20., ymin+3.*(ymax-ymin)/20., $
  'Cassini pointing to: ' + string(caspoint_lat, format='(d5.2)') + '!9%!XN, ' $
  + string(caspoint_lon, format='(d6.2)') + '!9%!XW', $
  color = red
;xyouts, xmin+(xmax-xmin)/20., ymin+2.*(ymax-ymin)/20., '3-db beam: 1.2!9%!X', color = red
xyouts, xmin+(xmax-xmin)/20., ymin+2.*(ymax-ymin)/20., 'HGA-FWHM: 2.4!9%!X', color = red
xyouts, xmin+(xmax-xmin)/20., ymin+(ymax-ymin)/20., 'HGA pointing error: !9+!X0.2!9%!X', color = red

xyouts, xmin+(xmax-xmin)/20., ymin+15.*(ymax-ymin)/20., wms, color = green

distgain = 20. * alog10(absc2h[0] / absc2h[2])
;xyouts, xmin+11.*(xmax-xmin)/20., ymin+17.*(ymax-ymin)/20., $
xyouts, xmin+13.*(xmax-xmin)/20., ymin+18.*(ymax-ymin)/20., $
  'distance gain since !18t!D0!N!X: ' + $
;  string(10. * alog10((absc2h[0] / absc2h[2])^2), format='(d4.1)') + ' dB', $
  string(distgain, format='(d4.1)') + ' dB', $
  color = orange

xyouts, xmin+11.*(xmax-xmin)/20., ymin+3.*(ymax-ymin)/20., $
  'range: ' + $
  string(absc2h[2], format = '(i5)') + ' km', $
  color = black

xyouts, xmin+11.*(xmax-xmin)/20., ymin+2.*(ymax-ymin)/20., $
  'range rate: ' + $
  string(-c2h_dot, format = '(d4.2)') + ' km/sec', $
  color = black

;xyouts, xmin+11.*(xmax-xmin)/20., ymin+18.*(ymax-ymin)/20., $
xyouts, xmin+11.*(xmax-xmin)/20., ymin+(ymax-ymin)/20., $
  'NCO frequency (mission mode): ' + $
  string(-c2h_dot / .299792458 * 2.040 - 38.51438, format = '(d5.1)') + ' kHz', $
;  string(-c2h_dot / .299792458 * 2.098 - 38.5, format = '(d5.1)') + ' kHz', $
  color = black

oplot, [xmin - 5, xmin - 5 - 5 * 2.225], [lat_ls0, lat_ls0], $
  color = black
xyouts, xmin - 5 - 2.5 * 2.225, lat_ls0 + .2, $
  '500 km', alignment = 0.5, color = black
oplot, [xmin - 5, xmin - 5], [lat_ls0 + 2.5 * 2.35, lat_ls0 - 2.5 * 2.35], $
  color = black
xyouts, xmin - 4.8, lat_ls0, $
  '500 km', alignment = 0.5, orientation = 90., color = black

;------------------ plot histogram of "good" missions ----------

;-----------------Include pointing error--------------
pe_flag = 1

if pe_flag then begin
  pe = 0.2 * randomn(seed, no_bases)  ; random number for HGA pointing error
  alpha = !dpi * randomu(seed, no_bases)  ; random number for HGA pointing error
  huy_offang = sqrt(huy_offang^2 + pe^2 - 2. * huy_offang * pe * cos(alpha))
  pe_string = '!C(included)'
endif else pe_string = '!C(not included)'
;-----------done: Include pointing error---------------

;huy_powdb[i] = 10 * log 2 * x^2 / x_3db^2    ; 3dB: 1.2 deg
huy_powdb = 3.0103 * huy_offang * huy_offang / 1.44
huy_powdb = 7.5 + distgain - huy_powdb

bins = 1
x = indgen(ceil((max(huy_powdb)-min(huy_powdb)) / bins)) * bins
x = x + floor(min(huy_powdb))
y = 100. * histogram(huy_powdb, binsize = bins) / no_bases

if n_elements(x) gt 0 and n_elements(y) gt 0 then begin
  xx = [x[0] - bins, x, x[n_elements(x) - 1] + bins]
  yy = [0., y, 0.]
endif else begin
  xx = [0., 1.]
  yy = [0., 1.]
endelse

plot, xx, yy, $
  psym = 10, $
  title = 'Link Margin Density Function at !18t-t!D0!N!X = ' + string(t_eval, format='(i3)') + ' min', $
  xtitle = 'SNR [dB]', $
  ytitle = 'PERCENTAGE (BINSIZE = ' + string(bins, format = '(i2)') + ' dB)', $
  xrange = [-20., 30.], $
  xstyle = 1, $
;  ystyle = 1, $
;  ystyle = 3, $
  color = black
oplot, [7.5, 7.5], [!y.crange[0], !y.crange[1]], linestyle = 1
oplot, [0., 0.], [!y.crange[0], !y.crange[1]], linestyle = 2

where_ok = where(xx ge 7.5, count)
if count ne 0 then $
  ok1 = total(yy[where_ok]) $
else ok1 = 0.0
where_ok = where(xx gt 0.0, count)
if count ne 0 then $
  ok2 = total(yy[where_ok]) $
else ok2 = 0.0
xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/25., $
        !y.crange[0]+11.*(!y.crange[1]-!y.crange[0])/20., $
  'Percentage of missions with!C1) link margin above 7.5 dB: ' + $
  string(ok1, format = '(d5.1)') + '%' + $
  '!C2) link margin above 0.0 dB: ' + $
  string(ok2, format = '(d5.1)') + '%'

xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/25., $
        !y.crange[0]+18.*(!y.crange[1]-!y.crange[0])/20., $
  'Link margin at !18t!D0!N!X: 7.5 dB', $
  color = black

xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/25., $
        !y.crange[0]+17.*(!y.crange[1]-!y.crange[0])/20., $
  'HGA pointing error: !9+!X0.2!9%!X' + pe_string, $
  color = black

xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/25., $
        !y.crange[0]+15.*(!y.crange[1]-!y.crange[0])/20., $
  'Huygens at !18t!D0!N!X: ' + string(lat_ls0, format='(d5.2)') + $
  '!9%!XN, ' + string(lon_ls0, format='(d6.2)') + '!9%!XW', $
  color = black

xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/25., $
        !y.crange[0]+14.*(!y.crange[1]-!y.crange[0])/20., $
  'Cassini pointing to: ' + string(caspoint_lat, format='(d5.2)') + '!9%!XN, ' $
  + string(caspoint_lon, format='(d6.2)') + '!9%!XW', $
  color = black

xyouts, !x.crange[0]+(!x.crange[1]-!x.crange[0])/25., $
        !y.crange[0]+13.*(!y.crange[1]-!y.crange[0])/20., $
  wms, $
  color = black

;------------------ done: plot histogram of "good" missions ----------

if (medium eq 0) or (medium eq 1) then device, /close
end