C+++++++++++++++++++++++++ C subroutine cappplc(alpha0, delta0, rmuafk, rmudfk, * parfk, rvfk, Jahr, Monat, itag, ihour, racioout, raequout, * decout) C C+++++++++++++++++++++++++ C computes the apparent place for a star at alpha0, delta0 C (J2000.0 degrees) with proper motions rmuafk, rmudfk (mas/yr, C cos(delta) applied to mu_alpha), with parallax parfk (arcsec) C and radial velocity rvfk, for year Jahr, month Monat, day itag, C and hour ihour. C The result will be returned in racioout (RA, CIO method), C raequout (RA, Equinox method) and decout (Declination for C both). implicit real*8 (a-h,o-z) dimension q(3),p1(3),p2(3),p3(3),xv(3),p4(3),p(3),PP(3) dimension EB(3),EV(3),SB(3),E(3),es(3) dimension r(6) real*8 m(3) real*8 RBPN(3,3),RC2I(3,3) character * 200 zeile C pi = 4.d0 * atan(1.d0) dtr = 180.d0/pi rmua0 = rmuafk / dtr / 3600.d03/ cos(delta0/dtr) * 100.d0 rmud0 = rmudfk / dtr / 3600.d03 * 100.d0 par = parfk / dtr / 3600.d00 C rv = rvfk * 21.095d0 13.02.12 rv = rvfk * 21.09495 d0 vaupar = par * rv C calculate Julian date call iau_CAL2JD(Jahr,Monat,itag,DJM0,DJ,ixflag) if(ixflag.ne.0) then write(6,777) ixflag 777 format(' Fehler bei CAL2JD Aufruf: ixflag = ',i6) stop end if DJM = DJM0 + DJ t1 = DJM + dble(ihour)/24.d0 t0 = 2451545.0 d0 C q(1) = cos(alpha0/dtr)*cos(delta0/dtr) q(2) = sin(alpha0/dtr)*cos(delta0/dtr) q(3) = sin(delta0/dtr) C write(6,'(f20.9,a)') q(1) , ' q' C write(6,'(f20.9,a)') q(2) , ' q' C write(6,'(f20.9,a)') q(3) , ' q' C m(1) = -cos(delta0/dtr)*sin(alpha0/dtr)*rmua0 & -sin(delta0/dtr)*cos(alpha0/dtr)*rmud0 & +cos(delta0/dtr)*cos(alpha0/dtr)*vaupar m(2) = cos(delta0/dtr)*cos(alpha0/dtr)*rmua0 & -sin(delta0/dtr)*sin(alpha0/dtr)*rmud0 & +cos(delta0/dtr)*sin(alpha0/dtr)*vaupar m(3) = cos(delta0/dtr) * rmud0 + sin(delta0/dtr)*vaupar C C write(6,'(f20.9,a)') C write(6,'(f20.9,a)') m(1) , ' m' C write(6,'(f20.9,a)') m(2) , ' m' C write(6,'(f20.9,a)') m(3) , ' m' C call DE405 barycentrischer positions- and velocity vector of earth: ET = T1 NTARG = 3 NCTR = 12 call pleph (ET, NTARG,NCTR,R) C write(6,'(f20.9,A)') R(1), ' R1-Erde' C write(6,'(f20.9,A)') R(2), ' R2-Erde' C write(6,'(f20.9,A)') R(3), ' R3-Erde' C write(6,'(f20.9,A)') R(4), ' V1-Erde' C write(6,'(f20.9,A)') R(5), ' V2-Erde' C write(6,'(f20.9,A)') R(6), ' V3-Erde' EB(1) = R(1) EB(2) = R(2) EB(3) = R(3) EV(1) = R(4) EV(2) = R(5) EV(3) = R(6) C C call DE405 barycentric position vector of sun: ET = T1 NTARG = 11 NCTR = 12 call pleph (ET, NTARG,NCTR,R) C write(6,'(f20.9,A)') R(1), ' R1-Sonne' C write(6,'(f20.9,A)') R(2), ' R2-Sonne' C write(6,'(f20.9,A)') R(3), ' R3-Sonne' SB(1) = R(1) SB(2) = R(2) SB(3) = R(3) C dt = (t1 - t0)/36525.d0 PP(1) = q(1) + dt*m(1) - par*EB(1) PP(2) = q(2) + dt*m(2) - par*EB(2) PP(3) = q(3) + dt*m(3) - par*EB(3) pscl = sqrt(PP(1)**2+PP(2)**2+PP(3)**2) p(1) = PP(1)/pscl p(2) = PP(2)/pscl p(3) = PP(3)/pscl C write(6,'(f20.9,a)') C write(6,'(f20.9,a)') p(1), ' p' C write(6,'(f20.9,a)') p(2), ' p' C write(6,'(f20.9,a)') p(3), ' p' C E(1) = EB(1) - SB(1) E(2) = EB(2) - SB(2) E(3) = EB(3) - SB(3) Escl = sqrt(E(1)**2 + E(2)**2 + E(3)**2) 01.08.13 escl ==> Escl es(1) = E(1)/Escl 01.08.13 es(2) = E(2)/Escl 01.08.13 es(3) = E(3)/Escl 01.08.13 C write(6,'(f20.9,a)') es(1), ' es(1)' C write(6,'(f20.9,a)') es(2), ' es(2)' C write(6,'(f20.9,a)') es(3), ' es(3)' C c Step 3 light deflection: C xmc2 = 9.87 D-09 16.02.12 C xmc2 = 9.8735 D-09 xmc2 = 9.8731 D-09 23.02.12 C call ldefl(es,p ,xmc2,p1) call ldefl(es,p ,xmc2,Escl,p1) 01.08.13 c C write(6,'(f20.9)') C write(6,'(f20.9,a)') p1(1), ' p1(1)' C write(6,'(f20.9,a)') p1(2), ' p1(2)' C write(6,'(f20.9,a)') p1(3), ' p1(3)' c Step 4 Aberration: call vframe(EV,p1,p2) c C write(6,'(f20.9)') C write(6,'(f20.9,a)') p2(1), ' p2(3)' C write(6,'(f20.9,a)') p2(2), ' p2(3)' C write(6,'(f20.9,a)') p2(3), ' p2(3)' c Step 5a (EQUINOX-Method) arg1 = t1 arg2 = 0.d0 CC call iau_PNM00A(arg1,arg2,RBPN) call iau_PNM06A(arg1,arg2,RBPN) 26.08.08 C write(6,'(f20.9,a)') RBPN(1,1),' 1-1 ' C write(6,'(f20.9,a)') RBPN(1,2),' 1-2 ' C write(6,'(f20.9,a)') RBPN(1,3),' 1-3 ' C p3(1) = RBPN(1,1)*p2(1)+RBPN(1,2)*p2(2)+RBPN(1,3)*p2(3) p3(2) = RBPN(2,1)*p2(1)+RBPN(2,2)*p2(2)+RBPN(2,3)*p2(3) p3(3) = RBPN(3,1)*p2(1)+RBPN(3,2)*p2(2)+RBPN(3,3)*p2(3) c C write(6,'(f20.9)') C write(6,'(f20.8,a)') p3(1), ' p3(1)' C write(6,'(f20.8,a)') p3(2), ' p3(2)' C write(6,'(f20.8,a)') p3(3), ' p3(3)' c Step 6a delta = asin(p3(3)) * dtr alpha = atan2(p3(2),p3(1)) * dtr if(alpha.lt.0d0) alpha = alpha + 360.d0 raequout = alpha C alpha = alpha / 15.d0 C idega = int(alpha) C imina= int(alpha*60.d0 - dble(idega)*60.d0) C seca = abs(alpha)*3600.d0 C & -dble(iabs(idega))*3600.d0-dble(imina)*60.d0 C write(6,999) idega , imina, seca 999 format(' alpha = ',I4,2x,i2,2x,f10.4,' Equinox-Method') decout = delta C idegd = int(delta) C imind = int(abs(delta*60.d0) - dble(iabs(idegd))*60.d0) C secd = abs(delta)*3600.d0 C & -dble(iabs(idegd))*3600.d0-dble(imind)*60.d0 C write(6,998) idegd , imind, secd 998 format(' delta = ',I4,2x,i2,2x,f10.4,' Equinox-Method') c Step 5b (CIO-Method) arg1 = t1 arg2 = 0.d0 C call iau_C2I00A(arg1,arg2,RC2I) call iau_C2I06A(arg1,arg2,RC2I) 26.08.08 C write(6,'(f20.9,a)') RC2I(1,1),' 1-1 ' C write(6,'(f20.9,a)') RC2I(1,2),' 1-2 ' C write(6,'(f20.9,a)') RC2I(1,3),' 1-3 ' C p3(1) = RC2I(1,1)*p2(1)+RC2I(1,2)*p2(2)+RC2I(1,3)*p2(3) p3(2) = RC2I(2,1)*p2(1)+RC2I(2,2)*p2(2)+RC2I(2,3)*p2(3) p3(3) = RC2I(3,1)*p2(1)+RC2I(3,2)*p2(2)+RC2I(3,3)*p2(3) c C write(6,'(f20.9)') C write(6,'(f20.9,a)') p3(1), ' p3(1)' C write(6,'(f20.9,a)') p3(2), ' p3(2)' C write(6,'(f20.9,a)') p3(3), ' p3(3)' c Step 6b delta = asin(p3(3)) * dtr alpha = atan2(p3(2),p3(1)) * dtr if(alpha.lt.0d0) alpha = alpha + 360.d0 racioout = alpha C alpha = alpha / 15.d0 C idega = int(alpha) C imina= int(alpha*60.d0 - dble(idega)*60.d0) C seca = abs(alpha)*3600.d0 C & -dble(iabs(idega))*3600.d0-dble(imina)*60.d0 C write(6,888) idega , imina, seca C imind = int(abs(delta*60.d0) - dble(iabs(idegd))*60.d0) C secd = abs(delta)*3600.d0 C & -dble(iabs(idegd))*3600.d0-dble(imind)*60.d0 C write(6,887) idegd , imind, secd 888 format(' alpha = ',I4,2x,i2,2x,f10.4,' CIO-Method') 887 format(' delta = ',I4,2x,i2,2x,f10.4,' CIO-Method') return end C subroutine ldefl(es,p1,xmc2,p2) subroutine ldefl(es,p1,xmc2,Escl,p2) 01.08.13 c c Computation of light deflection accoring to "Astronomical Almanac", S. B30 c --------------------------------------------------------------------------- c implicit real*8 (a-h,o-z) c dimension es(3),p1(3),p2(3) c c c Skalarprodukt bilden c call iau_pdp(p1,es,edotp) c c Betrag von ES bilden c esb = dsqrt(es(1)*es(1)+es(2)*es(2)+es(3)*es(3)) c c C p2(1) = p1(1) + 2.d0*xmc2*esb*(es(1)-edotp*p1(1))/(1.d0+edotp) C p2(2) = p1(2) + 2.d0*xmc2*esb*(es(2)-edotp*p1(2))/(1.d0+edotp) C p2(3) = p1(3) + 2.d0*xmc2*esb*(es(3)-edotp*p1(3))/(1.d0+edotp) p2(1) = p1(1) + 2.d0*xmc2/Escl*(es(1)-edotp*p1(1))/(1.d0+edotp) 01.08.13 p2(2) = p1(2) + 2.d0*xmc2/Escl*(es(2)-edotp*p1(2))/(1.d0+edotp) 01.08.13 p2(3) = p1(3) + 2.d0*xmc2/Escl*(es(3)-edotp*p1(3))/(1.d0+edotp) 01.08.13 c c return end subroutine vframe(xv,p2,p3) c c Transformation to the frame moving with velocity xv c ---------------------------------------------------- c implicit real*8 (a-h,o-z) C dimension xv(3),p2(3),p3(3) c C xvc = 0.0057755d0 16.02.12 C xvc = 0.005775518d0 xvc = 0.00577551833d0 23.02.12 c c c Transform V into light velocity. c V was expressed from the subroutine in au per day. c -------------------------------------------------- c c do i = 1,3 xv(i) = xvc*xv(i) end do c c c Berechnung des Betrages von V c xvb = dsqrt(xv(1)*xv(1)+xv(2)*xv(2)+xv(3)*xv(3)) c c C beta = 1.d0/dsqrt(1.d0 + (xv(1)*xv(1)+xv(2)*xv(2)+xv(3)*xv(3))) beta = 1.d0/dsqrt(1.d0 - (xv(1)*xv(1)+xv(2)*xv(2)+xv(3)*xv(3))) betam1 = 1.d0/beta c c Compute the scalar product c call iau_pdp (p2,xv,pdotv) c c x1 = 1.d0 + pdotv/(1.d0+betam1) do i = 1,3 C p3(i) = (betam1*p2(i)+ x1*xv(i))/(1+pdotv) 16.02.12 p3(i) = (betam1*p2(i)+ x1*xv(i))/(1.d0+pdotv) end do c c return end