WSJT-X/lib/azdist.f90

subroutine azdist(MyGrid,HisGrid,utch,nAz,nEl,nDmiles,nDkm,nHotAz,nHotABetter)

  character(len=*) :: MyGrid,HisGrid
  character*6 mygrid0,hisgrid0
  real*8 utch,utch0
  logical HotABetter,IamEast
  real eltab(22),daztab(22)
  data eltab/18.,15.,13.,11.,9.,8.,7.,6.,5.3,4.7,4.,3.3,2.7,    &
       2.,1.5,1.,0.8,0.6,0.4,0.2,0.0,0.0/
  data daztab/21.,18.,16.,15.,14.,13.,12.,11.,10.7,10.3,10.,    &
       10.,10.,10.,10.,10.,10.,9.,9.,9.,8.,8./
  data mygrid0/"      "/,hisgrid0/"      "/,utch0/-999.d0/
  save

  if(MyGrid.eq.HisGrid) then
     naz=0
     nel=0
     ndmiles=0
     ndkm=0
     nhotaz=0
     nhotabetter=1
     go to 999
  endif

  if(mygrid.eq.mygrid0 .and. hisgrid.eq.hisgrid0 .and.              &
       abs(utch-utch0).lt.0.1666667d0) go to 900
  utch0=utch
  mygrid0=mygrid
  hisgrid0=hisgrid
  utchours=utch

  if(MyGrid(5:5).eq.' ') MyGrid(5:5)='m'
  if(MyGrid(6:6).eq.' ') MyGrid(6:6)='m'
  if(HisGrid(5:5).eq.' ') HisGrid(5:5)='m'
  if(HisGrid(6:6).eq.' ') HisGrid(6:6)='m'

  if(MyGrid.eq.HisGrid) then
     Az=0.
     Dmiles=0.
     Dkm=0.0
     El=0.
     HotA=0.
     HotB=0.
     HotABetter=.true.
     go to 900
  endif
  call grid2deg(MyGrid,dlong1,dlat1)
  call grid2deg(HisGrid,dlong2,dlat2)
  eps=1.e-6
  Az=0.
  Dmiles=0.
  Dkm=0.0
  El=0.
  HotA=0.
  HotB=0.
  HotABetter=.true.
  if(abs(dlat1-dlat2).lt.eps .and. abs(dlong1-dlong2).lt.eps) go to 900

  difflong=mod(dlong1-dlong2+720.0,360.0)
  if(abs(dlat1+dlat2).lt.eps .and. abs(difflong-180.0).lt.eps) then
! Antipodes
     Dkm=20400
     go to 900
  endif

  call geodist(dlat1,dlong1,dlat2,dlong2,Az,Baz,Dkm)

  ndkm=Dkm/100
  j=ndkm-4
  if(j.lt.1) j=1
  if(j.gt.21)j=21
  if(Dkm.lt.500.0) then
     El=18.0
  else
     u=(Dkm-100.0*ndkm)/100.0
     El=(1.0-u)*eltab(j) + u*eltab(j+1)
  endif

  daz=daztab(j) + u * (daztab(j+1)-daztab(j))
  Dmiles=Dkm/1.609344

  tmid=mod(UTChours-0.5*(dlong1+dlong2)/15.0+48.0,24.0)
  IamEast=.false.
  if(dlong1.lt.dlong2) IamEast=.true.
  if(dlong1.eq.dlong2 .and. dlat1.gt.dlat2) IamEast=.false.
  azEast=baz
  if(IamEast) azEast=az
  if((azEast.ge.45.0 .and. azEast.lt.135.0) .or.                    &
       (azEast.ge.225.0 .and. azEast.lt.315.0)) then
! The path will be taken as "east-west".
     HotABetter=.true.
     if(abs(tmid-6.0).lt.6.0) HotABetter=.false.
     if((dlat1+dlat2)/2.0 .lt. 0.0) HotABetter=.not.HotABetter
  else
! The path will be taken as "north-south".
     HotABetter=.false.
     if(abs(tmid-12.0).lt.6.0) HotABetter=.true.
  endif
  if(IamEast) then
     HotA = Az - daz
     HotB = Az + daz
  else
     HotA = Az + daz
     HotB = Az - daz
  endif
  if(HotA.lt.0.0)   HotA=HotA+360.0
  if(HotA.gt.360.0) HotA=HotA-360.0
  if(HotB.lt.0.0)   HotB=HotB+360.0
  if(HotB.gt.360.0) HotB=HotB-360.0

900 continue
  naz=nint(Az)
  nel=nint(el)
  nDmiles=nint(Dmiles)
  nDkm=nint(Dkm)
  nHotAz=nint(HotB)
  nHotABetter=0
  if(HotABetter) then
     nHotAz=nint(HotA)
     nHotABetter=1
  endif

999 return
end subroutine azdist