subroutine orbitdelta (period, tper,ecc,axis,xincl,
     & xnode,xomega,
     & alpha0,delta0,rmuafk,rmudfk,parfk,rvfk,
     & jahr,monat,itag,ihour,dalf,ddel)
      implicit real*8 (a-h,o-z)
c
  
      pi = 4.d0 * atan(1.d0)
      trgr = pi/180.d0
c
      xnode_rad  = xnode  * trgr
      xomega_rad = xomega * trgr
      xincl_rad  = xincl  * trgr
      axis_sin   = axis / sin(xincl_rad)
c
      xao    = alpha0 * trgr
      xdo    = delta0 * trgr
      xmao   = rmuafk/3600.d03 * trgr / cos(xdo)
      xmdo   = rmudfk/3600.d03 * trgr
      paro  = parfk/3600.d0  * trgr                                             nur falls nach cappplc
      velo  =  rvfk
c     --------------------------------------------
      call iau_CAL2JD(Jahr,Monat,itag,DJM0F,DJF,ixflag)
      xjd = DJM0F + DJF + dble(ihour)/24.d0
      call jdje (xjd,time)                                                      !=> /progs/subprog.f   Jul.Datum => Jul. Epoche
C     write(*,*) 'time',time
      call binorb(period,tper,ecc,axis_sin,xincl_rad,
     f xnode_rad,xomega_rad,time,                 !=> Anhang
     f     xao,xdo,xmao,xmdo,paro,velo,dalf,ddel)
c     --------------------------------
c     ----------------------------------
c                                                                             !   in:   p4(3)    vector
c     Bahnkorrektion evtl. hinzurechnen                                       !   out:  theta    longitude [rad]
c                                                                             !         phi      latitude [rad]
C
      dalf  =   dalf/3600.d0         
      ddel  =   ddel/3600.d0           
C
c
      return                                                               !26.10.05    Ende Hauptprogramm
      end
      subroutine binorb(p,tp,e,ax,pinc,gom,om,t,
     f     xa,xd,xma,xmd,par,vel,xda,xdd)
c
      implicit real*8 (a-h,o-z)
c
      real*8 precmtx(3,3)
c
      pi  = 3.1415926535897932d0
      trg = pi/180.d0
c
c
      ln_out = 31
c
c
c
      ex = e
      pmu = 2.d0*pi/p
      pmea= pmu * (t - tp)
c     eps   = 1.0d-5
      eps   = 1.0d-10                                                   Wiel.-Jahr.29.11.99
      ge = datan2(dsin(pmea),dcos(pmea)-e)
      itn = 0
 10   ea = ge
      ge = pmea + e*dsin(ea)
      if(dabs(ge-ea).lt. eps) go to 20
      if(itn.gt. 80) then
c      pause ' itn > 80'
      go to 20
      end if
      go to 10
 20   radius = ax*(1.d0 - e * cos(ge))
      truano= 2.d0* atan(sqrt((1.d0+e)/(1.d0-e))* tan(ge*0.5d0))
      theta = atan2(sin(truano+om)*cos(pinc),cos(truano+om))+ gom
      rho   = radius * cos(truano+om)/cos(theta-gom)
      xda   = rho * sin(theta)
      xdd   = rho * cos(theta)
c
c     write(ln_out,'(//,''Korrektionen in mas  '',2f20.2)')  xda,xdd
c
c     Das sind die Korrektionen, bezogen auf das KoS J2000
c     Diese sind nun umzurechnen auf das Aequinox of Date.
c     Die Einheiten sind noch mas und in xda ist cos(delta) drin.
c
c     Umwandlung der Korrektionen in Radian,  cos(delta) eliminieren
c
      xda1 = xda
      xdd1 = xdd
c
      xda = (xda/3.6d6)*trg/dcos(xd)
      xdd = (xdd/3.6d6)*trg
c
c     write(ln_out,'(/,''Korrektionen in Radian  '',2f20.12)')
c    f   xda,xdd
c
c
c     Epoch of Date in Julianischem Datum berechnen
c     J2000  in Julianischem Datum berechnen
c
      call jejd (t,xt1)
      call jejd (2000.d0,xt0)
c     write(ln_out,'(//,''xt1 = '',2f20.3)')  xt0,xt1
c
c     Prazessionsmatrix rechnen
c
      call kprz(1,xt0,xt1,zet,teta,zeta)
      call przmtx(zet,teta,zeta,precmtx)
c
c     write(ln_out,'(''praez.matrix'',//,3(/,3f30.15))')  precmtx
c
c
c     Stern-Position (und Eigenbewegung) fuer Epoch of Date rechnen.
c     Aequinoktium bleibt.
c
c     write(ln_out,'(''Vor 1. Transformation '',2f20.3,/,4f20.10,/,
c    f    f20.10,f10.2)')  xt0,xt1,xa,xd,xma,xmd,par,vel
c     call tspvel (xt0,xt0,xt0,xt1,xa,xd,xma,xmd,par,vel,
c    f                             ya,yd,yma,ymd,ypar,yvel,ier)
c     cosd = dcos(yd)
c
c     write(ln_out,'(/,''Nach 1. Transformation '',/,2f20.10)') ya,yd
c
c     Nun die Doppelsternkorrektionen auf das Equinox of Date transformieren
c
c     write(ln_out,'(''Vor 2. Transformation '',2f20.3,/,4f20.10,/,
c    f    f20.10,f10.2)')  xt0,xt1,ya,yd,xda,xdd,ypar,yvel
c
      call preces(xa,xd,xda,xdd,precmtx,za,zd,zda,zdd)
c
c     call tspvel (xt0,xt1,xt1,xt1,ya,yd,xda,xdd,ypar,yvel,
c    f                             za,zd,zda,zdd,zpar,zvel,ier)
c     write(ln_out,'(/,''Nach 2. Transformation '',/,4f20.10)')
c    f          za,zd,zda,zdd
c
c     write(ln_out,'(''xda, xdd im bogenmass '',2f20.10)') zda,zdd
c
c     Ruecktransformation der umgerechneten Korrektionen in Bogensekunden,
c     wie sie vom APFS-Programm erwartet werden.
c
c     xda  = (zda/trg)*3600.d0/15.d0
      xda  = (zda/trg)*3600.d0
      xdd  = (zdd/trg)*3600.d0
c
c     Zum Vergleich nochmal alles in mas* rechnen
c
c     xda2 = xda*1000.d0*dcos(zd)*15.d0
c     xdd2 = xdd*1000.d0
c
c     xdiff = dsqrt(xda1*xda1+xdd1*xdd1) - dsqrt(xda2*xda2+xdd2*xdd2)
c
c     write(ln_out,'(f10.1,2f10.2,5x,3f10.2)') t,xda1,xda2,xdd1,xdd2,
c    f      xdiff
c
c
c     write(ln_out,'(//,
c    f   ''umgerechnete Korrektionen in Sec und arc-sec '',
c    f  2f20.8)')  xda,xdd
c
c
      return
      end