c .........................................................................
c Reduction to the HIPPARCOS system of positions and/or proper motions
c given in the system of one of the catalogues, discussed by Schwan
c in "Astronomy and Astrophysics, 373, 1099 (2001)" and in
c "Astronomy and Astrophysics, in press".
c
c
c The user has two options to use the program:
c
c 1) Full reduction to the HIPPARCOS system, including the global rotation
c and the elimination of regional errors.
c This option will in most cases meet the requirements of the user.
c
c 2) Reduction to the IAU (1976) system of Astronomical constants
c (no regional corrections or global rotations are applied).
c This option provides a catalogue which has the conventions of the
c FK5 catalogue, for example. It may be useful for those users who want to
c derive their own global rotation and/or regional corrections independently
c from those offered in the present software.
c
c
c
c What the user has to do:
c ========================
c
c 1) Provide the subroutine "liescat".
c This subroutine reads the catalogue.
c An example for such a subroutine with all necessary explanations
c is given in this program. Make the changes necessary to meet your
c catalogue format.
c
c Important Note: !!!!!!!!!
c If a magnitude equation exists in the catalogue comparison
c (i.e. there should be at least one function with p not equal 0
c for that comparison in the files koeff.pos or koeff.eb
c which define the comparison) then you should provide the
c HIPPARCOS magnitude as the star's apparent magnitude in the
c subroutine liescat for xmag, because the HIPPARCOS appaerent
c magnitude has been used in the determination of the
c systematic differences.
c 2) Arrange the output format in the subroutine "output"
c 3) The input catalogue should be in the file named "catorg.dat"
c the reduced catalogue will be in the file named "catred.dat"
c 4) All used file should be in the same directory. The files needed are:
c sydica.publ.f : FORTRAN Program. An executable file has to be
c produced by the user.
c catorg.dat : file with the input catalogue
c koeff.pos : file with the coefficients defining the rotations
c and systematic differences in positions.
c koeff.eb : file with the coefficients defining the rotations
c and systematic differences in proper motion.
c..............................................................................
c
c
implicit real*8 (a-h,o-z)
c
real*8 sysdif_eb(4),dsdelt_eb(4),dsalfa_eb(4),
f dsmag_eb(4),dsrot_eb(4),rstr_eb(4,2)
real*8 sysdif_pos(4),dsdelt_pos(4),dsalfa_pos(4),
f dsmag_pos(4),dsrot_pos(4),rstr_pos(4,2)
real*8 sysdif_colour_pos(4),dsdelt_colour_pos(4),
| dsalfa_colour_pos(4),dsmag_colour_pos(4),
| dsrot_colour_pos(4),rstr_colour_pos(4,2)
real*8 sysdif_colour_eb(4),dsdelt_colour_eb(4),
| dsalfa_colour_eb(4),dsmag_colour_eb(4),
| dsrot_colour_eb(4),rstr_colour_eb(4,2)
real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 sysreg(4)
c
character*12 a12,d12
c character feld1*500,feld2*160,feldh1*382
c
integer ivalid_pos(4),ivalid_eb(4)
integer ivalid_colour_pos(4),ivalid_colour_eb(4)
c
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb
character*1 answer,vorza,vorzd,minus,plus
character*50 catnames
c
dimension catnames(200)
c
c
common/sydpos/equi_pos,epa_pos,epd_pos,sysdif_pos,dsdelt_pos,
f dsalfa_pos,dsmag_pos,dsrot_pos,rstr_pos,
f ivalid_pos
common/sydeb /equi_eb ,epa_eb ,epd_eb ,sysdif_eb ,dsdelt_eb ,
f dsalfa_eb ,dsmag_eb ,dsrot_eb ,rstr_eb ,
f ivalid_eb
common/sydposcolour/equi_colour_pos,epa_colour_pos,epd_colour_pos,
f sysdif_colour_pos,dsdelt_colour_pos,
f dsalfa_colour_pos,dsmag_colour_pos,
f dsrot_colour_pos,rstr_colour_pos,ivalid_colour_pos
common/sydebcolour/equi_colour_eb,epa_colour_eb,epd_colour_eb,
f sysdif_colour_eb,dsdelt_colour_eb,
f dsalfa_colour_eb,dsmag_colour_eb,
f dsrot_colour_eb,rstr_colour_eb,ivalid_colour_eb
common/sigsys/sig_sysdif_pos,sig_sysdif_eb,
f sig_sysdif_colour_pos,sig_sysdif_colour_eb
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
common/xlimits/ xdmin,xdmax,xmmin,xmmax,xcmin,xcmax
c
pi = 4.d0*datan(1.d0)
tra = pi/12.d0
trd = pi/180.d0
c
c
plus = '+'
minus = '-'
c
c
20 continue
c
c Selection on the kind of reductions
c ===================================
c
c 1) full reduction to the HIPPARCOS system
c 2) Reduction merely to the IAU (1976) System of Astronomical Constants.
c
write(*,'(//,''Select the kind of reduction: '',//,
f ''1 = Full reduction to the Hipparcos system, including '',
f ''rotations and regional corrections'',/,
f ''2 = Reduction to the IAU (1976) System, without inclusion '',
f ''of rotations and regional corrections'',//,
f ''Please select '',/)')
read(*,'(i1)') icasesys
if (icasesys.eq.1) write(*,'(//,
f ''Full reduction to the Hipparcos system, including '',
f ''rotations and regional corrections'',//)')
if (icasesys.eq.2) write(*,'(//,
f ''Reduction to the IAU (1976) System, without inclusion '',
f ''of rotations and regional corrections'',//)')
write(*,'(//,''Is it o.k.? y/n'')')
read(*,'(a)') answer
if (answer.ne.'y') goto 20
c
c
25 continue
write(*,'(//,''Are colour dependent errors to be applied?'',//,
f ''1 = yes'',/,
f ''2 = no'')')
read(*,'(i1)') icasecol
if (icasecol.eq.1) write(*,'(/,
f ''Colour dependent errors will be applied'')')
if (icasecol.eq.2) write(*,'(/,
f ''Colour dependent errors will not be applied'')')
write(*,'(/,''Is that ok? y/n'')')
read(*,'(a)') answer
if (answer.ne.'y') goto 25
c
c
c Definition of Files:
c ====================
c
call open_files
c
c
c Definition of the system of the given observation
c =================================================
c
call select_syst(isyst)
c
c
c Definition of the catalogue names (for search in the coeeficient files);
c Definition of the conventions used to reduce to the HIPPARCOS system.
c ========================================================================
c
call define_names
c
call conventions(isyst,icsort,ip,equcat,eqk1,epeqk1,eqk2)
c
c
c Read the coefficients for the relevant catalogue comparison
c in the case icasesys = 1 (Full reduction to HIPPARCOS-system
c ==============================================================
c
write(*,'(''Files used: '',a,/,
f '' : '',a,/,
f '' : '',a,/,
f '' : '',a,//)')
f fnfk_pos,fnfk_eb,fnfk_pos,fnfk_eb
c
call sydica_pos(1,catnames(isyst),alpha,delta,xmag)
call sydica_eb (1,catnames(isyst),alpha,delta,xmag)
call sydica_colour_pos(1,catnames(isyst),alpha,delta,xmag)
call sydica_colour_eb (1,catnames(isyst),alpha,delta,xmag)
c
if (icasesys.eq.2) goto 65
c
xepf = 2000.0d0
xepv = epa_pos
c epverg = (epa_pos+epd_pos)/2.d0
deltat = (xepf-xepv)/100.d0
c
write(*,'(//,''Epoch for which the catalogue comparison '',
f ''holds:'',f10.3,/,
f ''Final epoch of the reduced catalogue '',
f '' :'',f10.3,//)')
f xepv,xepf
write(*,'(''Is it ok? y/n'',/)')
read(*,'(a)') answer
if (answer.eq.'y') goto 88
write(*,'(//,''Redefine epochs in the program '',//)')
goto 900
88 continue
c
65 continue
c
c
c .......................................
c Write heading line into the output file
c
c write(lnfk_red,'(
c f ''s-a = mean error of systematic difference'',
c f '' in right ascension'',/,
c f ''s-d = mean error of systematic difference'',
c f '' in declination '',/,
c f ''s-m-a = mean error of systematic difference'',
c f '' of proper motion in right ascension '',/,
c f ''s-m-d = mean error of systematic difference'',
c f '' of proper motion in declination '',/)')
c write(lnfk_red,'(
c f '' No. RA p.m. Dec. p.m. '',
c f '' mag s-a s-d s-m-a s-m-d'',/)')
c
c .........................................
c
rewind lnfk_orig
c
c
c read next star
c ==============
c
100 continue
c
call liescat(lnfk_orig,num,alpha,delta,xmya,xmyd,
f epa,epd,vel ,par,xmag,xcol,ier)
c
c
if (ier.lt.0) goto 900
c
c
c Check whether the star's declination and magnitude are within
c the limits of the catalogue comparison
c
call checklimits (num,delta,xmag,xcol,ier1,ier2,ier3)
if ( (ier1.eq.1) .and. (xmag.gt.xmmax) ) xmag = xmmax
if ( (ier1.eq.1) .and. (xmag.lt.xmmin) ) xmag = xmmin
if (ier3.eq.1) goto 100
c
c
c Computation of the systematic corrections for the star's
c right ascension alpha, declination delta, and apparent
c magnitude xmag;
c and eventually the colour equation.
c The determination of the global rotation and of the
c regional corrections has been made with a catalogue
c which was first reduced to the IAU (1976) system.
c This means in particular, that the coordinates used as
c arguments in those formulae were referred to the
c equinox and equator J2000. The epoch of the comparison
c is the mean epoch of the respective catalogue. All these
c have been read above from the coefficient-file.
c ========================================================
c
c
if (icasesys.eq.1) then
c
c
c Transformation of the position to the mean epoch and equinox
c to which the rotation parameters and the systematic differences
c are related.
c
c Transformation of epochs in Julian date
c
call jejd (equcat,eqcj)
call jejd (epa ,epcaj)
call jejd (epd ,epcdj)
call jejd (equi_pos,eqvposj)
call jejd (equi_eb ,eqvebj)
call jejd (epa_pos ,epvposj)
call jejd (epa_eb ,epvebj)
c
c Transformation of the position to the epoch and equinox for which the catlogue comparison holds
c
xa1 = alpha
xd1 = delta
xma1 = xmya
xmd1 = xmyd
xv1 = vel
xp1 = par
c
c
call tspv2(ip,eqcj,epcaj,eqvposj,epvposj,xa1,xd1 ,xma1, xmd1,
f xp1,xv1,
f xa2,xxx,xma2,yyyy,
f xp2,xv2,ier)
call tspv2(ip,eqcj,epcdj,eqvposj,epvposj,xa1,xd1 ,xma1, xmd1,
f xp1,xv1,
f xxx,xd2,yyyy,xmd2,
f xp2,xv2,ier)
c
c The positions for determining the systematic corrections in positions are:
c
xapos = xa2
xdpos = xd2
c
c
call tspv2(ip,eqcj,epcaj,eqvebj ,epvebj ,xa1,xd1 ,xma1, xmd1,
f xp1,xv1,
f xa2,xxx,xma2,yyyy,
f xp2,xv2,ier)
call tspv2(ip,eqcj,epcdj,eqvebj ,epvebj ,xa1,xd1 ,xma1, xmd1,
f xp1,xv1,
f xxx,xd2,yyyy,xmd2,
f xp2,xv2,ier)
c The positions for determining the systematic corrections in proper motion are:
c
xaeb = xa2
xdeb = xd2
c
c
c
call sydica_pos(2,catnames(isyst),xapos,xdpos,xmag)
call sydica_eb (2,catnames(isyst),xaeb ,xdeb ,xmag)
c
c
c Computation of colour dependent errors
c
sysdif_colour_pos(1) = 0.d0
sysdif_colour_pos(2) = 0.d0
sysdif_colour_eb (1) = 0.d0
sysdif_colour_eb (2) = 0.d0
if ( (icasecol.eq.1) .and. (ier3.eq.0) ) then
c
call sydica_colour_pos(2,catnames(isyst),xapos,xdpos,xcol)
call sydica_colour_eb (2,catnames(isyst),xaeb ,xdeb ,xcol)
c
c
end if
c
c
sysreg(1) = sysdif_pos(1) + sysdif_colour_pos(1)
sysreg(2) = sysdif_pos(2) + sysdif_colour_pos(2)
sysreg(3) = sysdif_eb (1) + sysdif_colour_eb(1)
sysreg(4) = sysdif_eb (2) + sysdif_colour_eb(2)
c
c
end if
c
c
c Define the rotations and regional corrections as zero in the
c case of casesys = 2, where merely the transition to the
c IAU (1976) System is made.
c ============================================================
c
if (icasesys.eq.2) then
do i = 1,4
sysreg(i) = 0.d0
c dsrot_pos(i) = 0.d0
c dsrot_eb(i) = 0.d0
end do
end if
c
c
call syst_trans (alpha,delta,xmya,xmyd,equcat,epa,epd,
f par,vel,ip,sysreg,xepv,
f eqk1,epeqk1,eqk2,icsort,
f ya,yd,yma,ymd)
c
alpha = ya
delta = yd
xmya = yma
xmyd = ymd
c
c Arrangement of the final output
c ===============================
c
call output (lnfk_red,num,alpha,delta,xmya,xmyd,epa,epd,xmag)
c
goto 100
c
900 stop
end
c
c
c
subroutine sydica_pos(icase,name_in,alpin,delin,amagin)
c
c Computation of systematic differences
c =====================================
c
c lnsyd: Unit reference number
c icase: 1 = reading and storing the coefficients at the beginning
c no systematic differences are computed.
c icase: 2 = computation of systematic differences at the position
c and for the magnitude of the star.
c isyst: defines the catalogue system and therefore the coefficients to be used
c alpin: Right ascension of the star in radians
c delin: Declination of the star in radians
c amag : Magnitude of the star in radians
c
c
c
c
c name Name of the catalogue (see subroutine define_catnames
c equi Equinox of the catalogue comparison.
c epa Epoch in R.A. for which the catalogue comparison holds
c epa Epoch in Dec. for which the catalogue comparison holds
c sysdif Total systematic differences in the same units as the coefficients
c dsdelt Part of systematic differences depending on the declination exlusively
c dsalfa Part of systematic differences depnding on Right ascension
c dsmag Part of systematic differences depnding on the Magnitude
c dsrot Global Rotation, eliminated in advance of the regional errors
c rstr dispersion of the the residulas and mean error of the system
c
C***********************************************************************
C
C
implicit real*8 (a-h,o-z)
C
character*50 name_in,name
c character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb
c
character feld*50
c
real*8 dsdelt_pos(4),dsalfa_pos(4),dsmag_pos(4),dsrot_pos(4),
f sysdif_pos(4),rstr_pos(4,2)
c real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 alpin,delin,pi,c1,c2,c3,c4
character*50 catnames
integer komp(4),jg(4),tk(4),h(50),jend(4),
f ivalid_pos(4),iexist(4)
dimension catnames(200)
c
real*8 danf(4),dend(4),norm(200,4),b(200,4),amaga(4),amage(4)
real*8 sb(200,4)
integer pj(200,4),nj(200,4),mj(200,4),lj(200,4),ind(200,4),p
c
real*8 normf,normh
c
common/sydpos/equi_pos,epa_pos,epd_pos,sysdif_pos,
f dsdelt_pos,dsalfa_pos,
f dsmag_pos,dsrot_pos,rstr_pos,ivalid_pos
c common/sigsys/sig_sysdif_pos,sig_sysdif_eb
c common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
c f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb
c common/xlimits/ xdmin,xdmax,xmmin,xmmax
c
common/sigsys/sig_sysdif_pos,sig_sysdif_eb,
f sig_sysdif_colour_pos,sig_sysdif_colour_eb
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
common/xlimits/ xdmin,xdmax,xmmin,xmmax,xcmin,xcmax
c
c
c
save
c
pi = 3.1415926535897932D0
c1=pi/180.d0
c2=pi/12.d0
c
c3=pi/(180.d0*3600.d0*100.d0)
c4=pi/(12.d0*3600.d0*1000.d0)
c
lnsyd = lnfk_pos
c
if( icase.ne. 1) goto 207
c
c Search of the required catalogue comparison and initialize the data
c
202 read(lnsyd,'(a)') name
if ( name_in.ne.name ) goto 202
c
read(lnsyd,'(4f8.2)') xdmin,xdmax,xmmin,xmmax
c
write(*,'(//,''The catalogue comparison can be applied for '',/,
f f8.2, '' < delta < '',f8.2)') xdmin,xdmax
write(*,'(/,''and for '',/,
f f8.2, '' < mag < '',f8.2)') xmmin,xmmax
c
read(lnsyd,'(i4,1x,3f9.3)',end=98)
* iak,equi_pos,epa_pos,epd_pos
c
c Read the rotational parameters.
c
read(lnsyd,'(3(2x,f8.2))') xom1,xom2,xom3
read(lnsyd,'(4f15.10)') a0,a1,b0,b1
c
c
do 50 i=1,iak
c
read(lnsyd,'(a)') feld
read(feld,'(i4)') l
komp(l) = l
read(feld,'(4x,2i6,2f14.5)')
* tk(l),jg(l),(rstr_pos(l,k),k=1,2)
c
jend(l) = jg(i)
je = jg(i)
c
c
c Computation of subscript defining the functions
c
do 60 j = 1,je
read(lnsyd,'(i10,f14.5,f15.4)') ind(j,l),b(j,l),sb(j,l)
call indux(ind(j,l),p,n,m,ll)
norm(j,l) = normh(p)*normf(n,m)
pj(j,l) = p
nj(j,l) = n
mj(j,l) = m
lj(j,l) = ll
60 continue
50 continue
c
return
c
207 continue
c
c end of catalogue search and data initialization
c
c Computation of systematic differences
c
xmagt = a0 + a1*amagin
alpha1 = alpin
delta1 = delin
sind = b0 + b1*dsin(delta1)
cosd = dcos(delta1)
C
do 100 k = 1,iak
C
sysdif_pos(K) = 0.d0
sig_sysdif_pos(k) = 0.d0
dsdelt_pos(K) = 0.d0
dsalfa_pos(K) = 0.d0
dsmag_pos (K) = 0.d0
dsrot_pos (K) = 0.d0
c
c Computation of infinitesimal rotation
c
call rotinf(alpha1,delta1,xom1,xom2,xom3,xarot,xdrot)
c
c
if (k.eq.1) dsrot_pos(k) = xarot*cosd
if (k.eq.2) dsrot_pos(k) = xdrot
c
c
c
do 120 j = 1,jend(k)
p = pj(j,k)
n = nj(j,k)
m = mj(j,k)
l = lj(j,k)
c
xsys = b(j,k)*norm(j,k)*xleg8(n,sind)
1 *four8(m,l,alpha1)*he(p,xmagt)
ysys = sb(j,k)*xsys/b(j,k)
c
c
sysdif_pos(k) = sysdif_pos(k) + xsys
sig_sysdif_pos(k) = sig_sysdif_pos(k) + ysys*ysys
if ( (m.eq.0) .and. (p.eq.0) )
f dsdelt_pos(k) = dsdelt_pos(k) + xsys
if ( (m.ne.0) .and. (p.eq.0) )
f dsalfa_pos(k) = dsalfa_pos(k) + xsys
if ( (p.ne.0) )
f dsmag_pos (k) = dsmag_pos(k) + xsys
120 continue
c
c Add rotation to the regional systematic error
c
sysdif_pos(k) = sysdif_pos(k) + dsrot_pos(k)
c
sig_sysdif_pos(k) = dsqrt(sig_sysdif_pos(k))
c
if (tk(k).ne.1) go to 125
125 continue
100 continue
c
return
c
c
c
98 continue
return
c
end
c
c
subroutine sydica_eb(icase,name_in,alpin,delin,amagin)
c
c Computation of systematic differences
c =====================================
c
c lnsyd: Unit reference number
c icase: 1 = reading and storing the coefficients at the beginning
c no systematic differences are computed.
c icase: 2 = computation of systematic differences at the position
c and for the magnitude of the star.
c isyst: defines the catalogue system and therefore the coefficients to be used
c alpin: Right ascension of the star in radians
c delin: Declination of the star in radians
c amag : Magnitude of the star in radians
c
c
c
c
c name Name of the catalogue (see subroutine define_catnames
c equi Equinox of the catalogue comparison.
c epa Epoch in R.A. for which the catalogue comparison holds
c epa Epoch in Dec. for which the catalogue comparison holds
c sysdif Total systematic differences in the same units as the coefficients
c dsdelt Part of systematic differences depending on the declination exlusively
c dsalfa Part of systematic differences depnding on Right ascension
c dsmag Part of systematic differences depnding on the Magnitude
c dsrot Global Rotation, eliminated in advance of the regional errors
c rstr dispersion of the the residulas and mean error of the system
c
c Requirements for a calling program:
c ===================================
c
c A calling program has to have the following statements:
c
c character*8 tel1,tel2
c real*8 dsdelt(4),dsalfa(4),dsmag(4),sysdif(4),rstr(4,2)
c real*8 alpin,delin
c integer ivalid(4)
c common/syd/equi,epa,epd,sysdif,dsdelt,dsalfa,dsmag,rstr,ivalid,
c f sigsys
c
c call sydica_eb(catnames,alpin,delin,amagin,lnsydk,lnsydl)
c
c
c***********************************************************************
c
c
implicit real*8 (a-h,o-z)
c
c
c character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb
character*50 name_in,name
c
character feld*50
c
c
real*8 dsdelt_eb(4),dsalfa_eb(4),dsmag_eb(4),dsrot_eb(4),
f sysdif_eb(4),rstr_eb(4,2)
c real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 alpin,delin,pi,c1,c2,c3,c4
character*50 catnames
dimension catnames(200)
integer komp(4),jg(4),tk(4),h(50),jend(4),
f ivalid_eb(4),iexist(4)
c
real*8 danf(4),dend(4),norm(200,4),b(200,4),amaga(4),amage(4)
real*8 sb(200,4)
c REAL*8 cov(4,200,200),cov1(200,200),fct(30000),sigsys(4)
integer pj(200,4),nj(200,4),mj(200,4),lj(200,4),ind(200,4),p
c
real*8 normf,normh
c
common/sydeb/equi_eb,epa_eb,epd_eb,sysdif_eb,dsdelt_eb,dsalfa_eb,
f dsmag_eb,dsrot_eb,rstr_eb,ivalid_eb
c common/sigsys/sig_sysdif_pos,sig_sysdif_eb
common/sigsys/sig_sysdif_pos,sig_sysdif_eb,
f sig_sysdif_colour_pos,sig_sysdif_colour_eb
c common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
c f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
c common/xlimits/ xdmin,xdmax,xmmin,xmmax
common/xlimits/ xdmin,xdmax,xmmin,xmmax,xcmin,xcmax
c
c
save
c
pi = 3.1415926535897932D0
c1=pi/180.d0
c2=pi/12.d0
c
c3=pi/(180.d0*3600.d0*100.d0)
c4=pi/(12.d0*3600.d0*1000.d0)
c
lnsyd = lnfk_eb
C
C
if( icase.ne. 1) goto 207
c
c Search of the required catalogue comparison and initialize the data
C
c
rewind lnsyd
c
202 read(lnsyd,'(a)') name
c
if ( name_in.ne.name ) goto 202
c
read(lnsyd,'(4f8.2)') xdmin,xdmax,xmmin,xmmax
c
write(*,'(//,''The catalogue comparison can be applied for '',/,
f f8.2, '' < delta < '',f8.2)') xdmin,xdmax
write(*,'(/,''and for '',/,
f f8.2, '' < mag < '',f8.2)') xmmin,xmmax
c
c
c The catalogue required has been found
c Read input data
c
c
read(lnsyd,'(i4,1x,3f9.3)',end=98)
* iak,equi_eb,epa_eb,epd_eb
c
c Read the rotational parameters.
c
read(lnsyd,'(3(2x,f8.2))') xom1,xom2,xom3
read(lnsyd,'(4f15.10)') a0,a1,b0,b1
c
do 50 i=1,iak
c
read(lnsyd,'(a)') feld
read(feld,'(i4)') l
komp(l) = l
read(feld,'(4x,2i6,2f14.5)')
* tk(l),jg(l),(rstr_eb(l,k),k=1,2)
c
jend(l) = jg(i)
je = jg(i)
c
c Computation of subscripts defining the functions
c
do 60 j = 1,je
read(lnsyd,'(i10,f14.5,f15.3)') ind(j,l),b(j,l),sb(j,l)
call indux(ind(j,l),p,n,m,ll)
norm(j,l) = normh(p)*normf(n,m)
pj(j,l) = p
nj(j,l) = n
mj(j,l) = m
lj(j,l) = ll
60 continue
50 continue
c
return
c
207 continue
c
c end of catalogue search and data initialization
c
c Computation of systematic differences
c
c
xmagt = a0 + a1*amagin
alpha1 = alpin
delta1 = delin
sind = b0 + b1*dsin(delta1)
cosd = dcos(delta1)
C
do 100 k = 1,iak
C
c
sysdif_eb(K) = 0.d0
sig_sysdif_eb(K) = 0.d0
dsdelt_eb(K) = 0.d0
dsalfa_eb(K) = 0.d0
dsmag_eb (K) = 0.d0
dsrot_eb (K) = 0.d0
c
c Computation of infinitesimal rotation
c
call rotinf(alpha1,delta1,xom1,xom2,xom3,xarot,xdrot)
if (k.eq.1) dsrot_eb(k) = xarot*cosd
if (k.eq.2) dsrot_eb(k) = xdrot
c
c
do 120 j = 1,jend(k)
p = pj(j,k)
n = nj(j,k)
m = mj(j,k)
l = lj(j,k)
c
xsys = b(j,k)*norm(j,k)*xleg8(n,sind)
1 *four8(m,l,alpha1)*he(p,xmagt)
ysys = sb(j,k)*xsys/b(j,k)
c
c
sysdif_eb(k) = sysdif_eb(k) + xsys
sig_sysdif_eb(k) = sig_sysdif_eb(k) + ysys*ysys
c
if ( (m.eq.0) .and. (p.eq.0) )
f dsdelt_eb(k) = dsdelt_eb(k) + xsys
if ( (m.ne.0) .and. (p.eq.0) )
f dsalfa_eb(k) = dsalfa_eb(k) + xsys
if ( (p.ne.0) )
f dsmag_eb (k) = dsmag_eb (k) + xsys
120 continue
c
c Add rotation to the regional systematic errors
c
sysdif_eb(k) = sysdif_eb(k) + dsrot_eb(k)
c
sig_sysdif_eb(k) = dsqrt(sig_sysdif_eb(k))
c
if (tk(k).ne.1) go to 125
125 continue
c end if
100 continue
c
return
c
c
c
98 continue
return
c
return
c
end
c
c
subroutine sydica_colour_pos(icase,name_in,alpin,delin,amagin)
c
c Computation of systematic differences
c =====================================
c
c lnsyd: Unit reference number
c icase: 1 = reading and storing the coefficients at the beginning
c no systematic differences are computed.
c icase: 2 = computation of systematic differences at the position
c and for the magnitude of the star.
c isyst: defines the catalogue system and therefore the coefficients to be used
c alpin: Right ascension of the star in radians
c delin: Declination of the star in radians
c amag : Colour index of the star
c
c
c
c
c name Name of the catalogue (see subroutine define_catnames
c equi Equinox of the catalogue comparison.
c epa Epoch in R.A. for which the catalogue comparison holds
c epa Epoch in Dec. for which the catalogue comparison holds
c sysdif Total systematic differences in the same units as the coefficients
c dsdelt Part of systematic differences depending on the declination exlusively
c dsalfa Part of systematic differences depnding on Right ascension
c dsmag Part of systematic differences depnding on the Magnitude
c dsrot Global Rotation, eliminated in advance of the regional errors
c rstr dispersion of the the residulas and mean error of the system
c
C***********************************************************************
C
C
implicit real*8 (a-h,o-z)
C
character*50 name_in,name
c character*50 fnfk_orig,fnfk_red,fnfk_colour_pos,fnfk_colour_eb
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb
c
character feld*50
c
real*8 dsdelt_colour_pos(4),dsalfa_colour_pos(4),
f dsmag_colour_pos(4),dsrot_colour_pos(4),
f sysdif_colour_pos(4),rstr_colour_pos(4,2)
c real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 alpin,delin,pi,c1,c2,c3,c4
character*50 catnames
integer komp(4),jg(4),tk(4),h(50),jend(4),
f ivalid_colour_pos(4),iexist(4)
dimension catnames(200)
c
real*8 danf(4),dend(4),norm(200,4),b(200,4),amaga(4),amage(4)
real*8 sb(200,4)
integer pj(200,4),nj(200,4),mj(200,4),lj(200,4),ind(200,4),p
c
real*8 normf,normh
c
c common/sydpos/equi_pos,epa_pos,epd_pos,sysdif_pos,
c f dsdelt_pos,dsalfa_pos,
c f dsmag_pos,dsrot_pos,rstr_pos,ivalid_pos
c common/sigsys/sig_sysdif_pos,sig_sysdif_eb
c common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
c f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb
c common/xlimits/ xdmin,xdmax,xmmin,xmmax
c
common/sydposcolour/equi_colour_pos,epa_colour_pos,epd_colour_pos,
f sysdif_colour_pos,dsdelt_colour_pos,
f dsalfa_colour_pos,dsmag_colour_pos,
f dsrot_colour_pos,rstr_colour_pos,ivalid_colour_pos
common/sigsys/sig_sysdif_pos,sig_sysdif_eb,
f sig_sysdif_colour_pos,sig_sysdif_colour_eb
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
common/xlimits/ xdmin,xdmax,xmmin,xmmax,xcmin,xcmax
c
c
save
c
pi = 3.1415926535897932D0
c1=pi/180.d0
c2=pi/12.d0
c
c3=pi/(180.d0*3600.d0*100.d0)
c4=pi/(12.d0*3600.d0*1000.d0)
c
lnsyd = lnfk_colour_pos
c
if( icase.ne. 1) goto 207
c
c Search of the required catalogue comparison and initialize the data
c
202 read(lnsyd,'(a)') name
if ( name_in.ne.name ) goto 202
c
read(lnsyd,'(4f8.2)') xdmin,xdmax,xmmin,xmmax
c
write(*,'(//,''The catalogue comparison can be applied for '',/,
f f8.2, '' < delta < '',f8.2)') xdmin,xdmax
write(*,'(/,''and for '',/,
f f8.2, '' < mag < '',f8.2)') xmmin,xmmax
c
read(lnsyd,'(i4,1x,3f9.3)',end=98)
* iak,equi_colour_pos,epa_colour_pos,epd_colour_pos
c
c Read the rotational parameters.
c
read(lnsyd,'(3(2x,f8.2))') xom1,xom2,xom3
read(lnsyd,'(4f15.10)') a0,a1,b0,b1
c
c
do 50 i=1,iak
c
read(lnsyd,'(a)') feld
read(feld,'(i4)') l
komp(l) = l
read(feld,'(4x,2i6,2f14.5)')
* tk(l),jg(l),(rstr_colour_pos(l,k),k=1,2)
ccc write(* ,'(4x,2i6,2f14.5)')
ccc * tk(l),jg(l),(rstr_colour_pos(l,k),k=1,2)
c
jend(l) = jg(i)
je = jg(i)
c
c
c Computation of subscript defining the functions
c
do 60 j = 1,je
read(lnsyd,'(i10,f14.5,f15.4)') ind(j,l),b(j,l),sb(j,l)
call indux(ind(j,l),p,n,m,ll)
norm(j,l) = normh(p)*normf(n,m)
pj(j,l) = p
nj(j,l) = n
mj(j,l) = m
lj(j,l) = ll
60 continue
50 continue
c
return
c
207 continue
c
c end of catalogue search and data initialization
c
c Computation of systematic differences
c
xmagt = a0 + a1*amagin
alpha1 = alpin
delta1 = delin
sind = b0 + b1*dsin(delta1)
cosd = dcos(delta1)
C
do 100 k = 1,iak
C
sysdif_colour_pos(K) = 0.d0
sig_sysdif_colour_pos(k) = 0.d0
dsdelt_colour_pos(K) = 0.d0
dsalfa_colour_pos(K) = 0.d0
dsmag_colour_pos (K) = 0.d0
dsrot_colour_pos (K) = 0.d0
c
c Computation of infinitesimal rotation
c
c call rotinf(alpha1,delta1,xom1,xom2,xom3,xarot,xdrot)
c
c
c if (k.eq.1) dsrot_pos(k) = xarot*cosd
c if (k.eq.2) dsrot_pos(k) = xdrot
c
c
c
do 120 j = 1,jend(k)
p = pj(j,k)
n = nj(j,k)
m = mj(j,k)
l = lj(j,k)
c
xsys = b(j,k)*norm(j,k)*xleg8(n,sind)
1 *four8(m,l,alpha1)*he(p,xmagt)
ysys = sb(j,k)*xsys/b(j,k)
c
c
sysdif_colour_pos(k) = sysdif_colour_pos(k) + xsys
sig_sysdif_colour_pos(k) = sig_sysdif_colour_pos(k) + ysys*ysys
if ( (m.eq.0) .and. (p.eq.0) )
f dsdelt_colour_pos(k) = dsdelt_colour_pos(k)+xsys
if ( (m.ne.0) .and. (p.eq.0) )
f dsalfa_colour_pos(k) = dsalfa_colour_pos(k)+xsys
if ( (p.ne.0) )
f dsmag_colour_pos (k) = dsmag_colour_pos(k)+xsys
120 continue
c
c Add rotation to the regional systematic error
c
sysdif_colour_pos(k) = sysdif_colour_pos(k) + dsrot_colour_pos(k)
c
sig_sysdif_colour_pos(k) = dsqrt(sig_sysdif_colour_pos(k))
c
if (tk(k).ne.1) go to 125
125 continue
100 continue
c
return
c
c
c
98 continue
return
c
end
c
c
subroutine sydica_colour_eb(icase,name_in,alpin,delin,amagin)
c
c Computation of systematic differences
c =====================================
c
c lnsyd: Unit reference number
c icase: 1 = reading and storing the coefficients at the beginning
c no systematic differences are computed.
c icase: 2 = computation of systematic differences at the position
c and for the magnitude of the star.
c isyst: defines the catalogue system and therefore the coefficients to be used
c alpin: Right ascension of the star in radians
c delin: Declination of the star in radians
c amag : Colour index of the star
c
c
c
c
c name Name of the catalogue (see subroutine define_catnames
c equi Equinox of the catalogue comparison.
c epa Epoch in R.A. for which the catalogue comparison holds
c epa Epoch in Dec. for which the catalogue comparison holds
c sysdif Total systematic differences in the same units as the coefficients
c dsdelt Part of systematic differences depending on the declination exlusively
c dsalfa Part of systematic differences depnding on Right ascension
c dsmag Part of systematic differences depnding on the Magnitude
c dsrot Global Rotation, eliminated in advance of the regional errors
c rstr dispersion of the the residulas and mean error of the system
c
C***********************************************************************
C
C
implicit real*8 (a-h,o-z)
C
character*50 name_in,name
c character*50 fnfk_orig,fnfk_red,fnfk_colour_pos,fnfk_colour_eb
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb
c
character feld*50
c
real*8 dsdelt_colour_eb(4),dsalfa_colour_eb(4),
f dsmag_colour_eb(4),dsrot_colour_eb(4),
f sysdif_colour_eb(4),rstr_colour_eb(4,2)
c real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
real*8 sig_sysdif_colour_pos(4),sig_sysdif_colour_eb(4)
real*8 alpin,delin,pi,c1,c2,c3,c4
character*50 catnames
integer komp(4),jg(4),tk(4),h(50),jend(4),
f ivalid_colour_eb(4),iexist(4)
dimension catnames(200)
c
real*8 danf(4),dend(4),norm(200,4),b(200,4),amaga(4),amage(4)
real*8 sb(200,4)
integer pj(200,4),nj(200,4),mj(200,4),lj(200,4),ind(200,4),p
c
real*8 normf,normh
c
c common/sydpos/equi_pos,epa_pos,epd_pos,sysdif_pos,
c f dsdelt_pos,dsalfa_pos,
c f dsmag_pos,dsrot_pos,rstr_pos,ivalid_pos
c common/sigsys/sig_sysdif_pos,sig_sysdif_eb
c common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
c f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb
c common/xlimits/ xdmin,xdmax,xmmin,xmmax
c
common/sydebcolour/equi_colour_eb,epa_colour_eb,epd_colour_eb,
f sysdif_colour_eb,dsdelt_colour_eb,
f dsalfa_colour_eb,dsmag_colour_eb,
f dsrot_colour_eb,rstr_colour_eb,ivalid_colour_eb
common/sigsys/sig_sysdif_pos,sig_sysdif_eb,
f sig_sysdif_colour_pos,sig_sysdif_colour_eb
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
common/xlimits/ xdmin,xdmax,xmmin,xmmax,xcmin,xcmax
c
c
save
c
pi = 3.1415926535897932D0
c1=pi/180.d0
c2=pi/12.d0
c
c3=pi/(180.d0*3600.d0*100.d0)
c4=pi/(12.d0*3600.d0*1000.d0)
c
lnsyd = lnfk_colour_eb
c
if( icase.ne. 1) goto 207
c
c Search of the required catalogue comparison and initialize the data
c
202 read(lnsyd,'(a)') name
if ( name_in.ne.name ) goto 202
c
read(lnsyd,'(4f8.2)') xdmin,xdmax,xmmin,xmmax
c
write(*,'(//,''The catalogue comparison can be applied for '',/,
f f8.2, '' < delta < '',f8.2)') xdmin,xdmax
write(*,'(/,''and for '',/,
f f8.2, '' < mag < '',f8.2)') xmmin,xmmax
c
read(lnsyd,'(i4,1x,3f9.3)',end=98)
* iak,equi_colour_eb,epa_colour_eb,epd_colour_eb
c
c Read the rotational parameters.
c
read(lnsyd,'(3(2x,f8.2))') xom1,xom2,xom3
read(lnsyd,'(4f15.10)') a0,a1,b0,b1
c
c
do 50 i=1,iak
c
read(lnsyd,'(a)') feld
read(feld,'(i4)') l
komp(l) = l
read(feld,'(4x,2i6,2f14.5)')
* tk(l),jg(l),(rstr_colour_eb(l,k),k=1,2)
c
jend(l) = jg(i)
je = jg(i)
c
c
c Computation of subscript defining the functions
c
do 60 j = 1,je
read(lnsyd,'(i10,f14.5,f15.4)') ind(j,l),b(j,l),sb(j,l)
call indux(ind(j,l),p,n,m,ll)
norm(j,l) = normh(p)*normf(n,m)
pj(j,l) = p
nj(j,l) = n
mj(j,l) = m
lj(j,l) = ll
60 continue
50 continue
c
return
c
207 continue
c
c end of catalogue search and data initialization
c
c Computation of systematic differences
c
xmagt = a0 + a1*amagin
alpha1 = alpin
delta1 = delin
sind = b0 + b1*dsin(delta1)
cosd = dcos(delta1)
C
do 100 k = 1,iak
C
sysdif_colour_eb(K) = 0.d0
sig_sysdif_colour_eb(k) = 0.d0
dsdelt_colour_eb(K) = 0.d0
dsalfa_colour_eb(K) = 0.d0
dsmag_colour_eb (K) = 0.d0
dsrot_colour_eb (K) = 0.d0
c
c Computation of infinitesimal rotation
c
c call rotinf(alpha1,delta1,xom1,xom2,xom3,xarot,xdrot)
c
c
c if (k.eq.1) dsrot_pos(k) = xarot*cosd
c if (k.eq.2) dsrot_pos(k) = xdrot
c
c
c
do 120 j = 1,jend(k)
p = pj(j,k)
n = nj(j,k)
m = mj(j,k)
l = lj(j,k)
c
xsys = b(j,k)*norm(j,k)*xleg8(n,sind)
1 *four8(m,l,alpha1)*he(p,xmagt)
ysys = sb(j,k)*xsys/b(j,k)
c
c
sysdif_colour_eb(k) = sysdif_colour_eb(k) + xsys
sig_sysdif_colour_eb(k) = sig_sysdif_colour_eb(k) + ysys*ysys
if ( (m.eq.0) .and. (p.eq.0) )
f dsdelt_colour_eb(k) = dsdelt_colour_eb(k)+xsys
if ( (m.ne.0) .and. (p.eq.0) )
f dsalfa_colour_eb(k) = dsalfa_colour_eb(k)+xsys
if ( (p.ne.0) )
f dsmag_colour_eb (k) = dsmag_colour_eb(k)+xsys
120 continue
c
c Add rotation to the regional systematic error
c
sysdif_colour_eb(k) = sysdif_colour_eb(k) + dsrot_colour_eb(k)
c
sig_sysdif_colour_eb(k) = dsqrt(sig_sysdif_colour_eb(k))
c
if (tk(k).ne.1) go to 125
125 continue
100 continue
c
return
c
c
c
98 continue
return
c
end
c
c
c
subroutine select_syst(isyst)
c
c Definition of the system of the given observation
c
implicit real*8 (a-h,o-z)
c
character answer
c
100 continue
c
write(*,'(''In which system is the observation given?'',//,
f '' 1 = FK5 Basic '',/,
f '' 2 = FK5 (Bright Extension + Sup)'',/,
f '' 3 = FK5 (Faint Extension) '',/,
f '' 4 = FK4 '',/,
f '' 5 = FK3 '',/,
f '' 6 = NFK '',/,
f '' 7 = N30 '',/,
f '' 8 = GC '',/,
f '' 9 = PGC '',/,
f ''10 = PPM '',/,
f ''11 = Eichelberger '',/,
f ''12 = SAO '',/,
f ''13 = AGK3 '',/,
f ''14 = AGK3R'',/,
f ''15 = Perth 70'',/,
f ''16 = IRS (FK4-based version)'',/,
f ''17 = IRS (FK5-based version)'',/,
f ''18 = AGK2 printed version'',/,
f ''19 = AGK2 revised version'',/,
f ''20 = AGK2A'',/,
f ''21 = ACRS (FK5-based)'',/,
f ''22 = SRS (FK4 based version)'',/,
f ''23 = SRS (FK5 based version)'',/,
f ''24 = CPC2 (FK5 based version)'',/,
f ''Please select '',/)')
read(*,'(i2)') isyst
if (isyst.eq.1)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the FK5 Basic system'',/)')
if (isyst.eq.2)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the FK5 (Bright Extension + Sup) system'',/)')
if (isyst.eq.3)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the FK5 (Faint Extension) system'',/)')
if (isyst.eq.4)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the FK4 system'',/)')
if (isyst.eq.5)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the FK3 system'',/)')
if (isyst.eq.6)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the NFK system'',/)')
if (isyst.eq.7)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the N30 system'',/)')
if (isyst.eq.8)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the GC system'',/)')
if (isyst.eq.9)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the PGC system'',/)')
if (isyst.eq.10)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the PPM system'',/)')
if (isyst.eq.11)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of Eichelbergers catalogue'',/)')
if (isyst.eq.12)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the SAO system'',/)')
if (isyst.eq.13)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the AGK3 system'',/)')
if (isyst.eq.14)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the AGK3R system'',/)')
if (isyst.eq.15)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the Perth 70 catalogue'',/)')
if (isyst.eq.16)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the IRS (FK4-based version) '',/)')
if (isyst.eq.17)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the IRS (FK5-based version) '',/)')
if (isyst.eq.18)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the AGK2 (printed version) '',/)')
if (isyst.eq.19)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the AGK2 (revised version) '',/)')
if (isyst.eq.20)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the AGK2A '',/)')
if (isyst.eq.21)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the ACRS(on FK5) '',/)')
if (isyst.eq.22)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the SRS (FK4-based version)'',/)')
if (isyst.eq.23)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the SRS (FK5-based version)'',/)')
if (isyst.eq.24)
f write(*,'(//,''It is assumed that the position is given '',
f ''in the system of the CPC2 (FK5-based version)'',/)')
c
write(*,'(''Is it o.k. y/n'',/)')
read(*,'(a)') answer
if (answer.ne.'y') goto 100
c
return
end
c
subroutine open_files
c
implicit real*8 (a-h,o-z)
c
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb
character*50 fnfk_colour_pos,fnfk_colour_eb
character*1 answer
character*50 catnames
dimension catnames(200)
c
c common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
c f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb
c
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
c
c
c Definition of used files:
c
c fk_orig = File with the given observation (Input file)
c fk_red = File with the reduce observation (Output file)
c fnfk_pos = File with the coefficients of regional differences in position
c fnfk_eb = File with the coefficients of regional differences in prop.mot.
c fnfk_colour_pos = File with the coefficients of the colour equatiposition
c in position
c fnfk_colour_eb = File with the coefficients of the colour equatiposition
c in declination
c
c
lnfk_orig = 11
lnfk_red = 12
lnfk_pos = 13
lnfk_eb = 14
lnfk_colour_pos = 15
lnfk_colour_eb = 16
c
c
fnfk_orig = 'catorg.dat'
fnfk_red = 'catred.dat'
fnfk_pos = 'koeff.pos'
fnfk_eb = 'koeff.eb'
fnfk_colour_pos = 'koeff.colour.pos'
fnfk_colour_eb = 'koeff.colour.eb'
c
300 continue
c
write(*,'(//,''The following files are used:'',//,
f ''Original observations : '',a,/,
f ''Transformed Observations: : '',a,/,
f ''Coefficients of systematic differences in position: : '',a,/,
f ''Coefficients of systematic differences in proper motion: '',a,/,
f ''Coefficients of colour equation in position : '',a,/,
f ''Coefficients of colour equation in proper motion : '',a,
f //,''Is it o.k.? y/n'' ,/)')
f fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb
read(*,'(a)') answer
if (answer.eq.'y') goto 400
c
write(*,'(/,''Definition of used files '',/)')
c
write(*,'(/,''Write file-name with original observations'')')
read(*,'(a)') fnfk_orig
write(*,'(/,''Write file-name with reduced observations'')')
read(*,'(a)') fnfk_red
write(*,'(/,''Write file-name with coeff. for positions '')')
read(*,'(a)') fnfk_pos
write(*,'(/,''Write file-name with coeff. for prop.mot. '')')
read(*,'(a)') fnfk_eb
c
goto 300
c
400 continue
c
c Opening of the files
c
open (lnfk_orig ,file=fnfk_orig ,status='old')
open (lnfk_red ,file=fnfk_red )
open (lnfk_pos ,file=fnfk_pos ,status='old')
open (lnfk_eb ,file=fnfk_eb ,status='old')
open (lnfk_colour_pos ,file=fnfk_colour_pos ,status='old')
open (lnfk_colour_eb ,file=fnfk_colour_eb ,status='old')
c
return
end
c
C
subroutine define_names
c
implicit real*8 (a-h,o-z)
c
character*50 catnames
character*50 fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb
character*50 fnfk_colour_pos,fnfk_colour_eb
c
dimension catnames(200)
c
c common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
c f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb
c
common/both/catnames,fnfk_orig,fnfk_red,fnfk_pos,fnfk_eb,
f lnfk_orig,lnfk_red,lnfk_pos,lnfk_eb,
f fnfk_colour_pos,fnfk_colour_eb,
f lnfk_colour_pos,lnfk_colour_eb
c
catnames( 1) = 'FK5 (Basic)'
catnames( 2) = 'FK5 (Bright Extension + Sup)'
catnames( 3) = 'FK5 (Faint Extension)'
catnames( 4) = 'FK4'
catnames( 5) = 'FK3'
catnames( 6) = 'NFK'
catnames( 7) = 'N30'
catnames( 8) = 'GC'
catnames( 9) = 'PGC'
catnames(10) = 'PPM'
catnames(11) = 'Eichelberger'
catnames(12) = 'SAO'
catnames(13) = 'AGK3'
catnames(14) = 'AGK3R'
catnames(15) = 'Perth 70 (Total catalogue)'
catnames(16) = 'IRS (FK4-based version)'
catnames(17) = 'IRS (FK5-based version)'
catnames(18) = 'AGK2 printed version'
catnames(19) = 'AGK2 revised version'
catnames(20) = 'AGK2A'
catnames(21) = 'ACRS (FK5-based version)'
catnames(22) = 'SRS (FK4-based version)'
catnames(23) = 'SRS (FK5-based version)'
catnames(24) = 'CPC2 (FK5-based version)'
c
return
end
c
subroutine conventions(isyst,icsort,ip,equcat,eqk1,epeqk1,eqk2)
c
c Definition of global conventions for the Input-Catalogue
c These are the following parameters which are the same for
c all stars in the catalogue.
c
c isyst : Number defining the catalgue system (chosen in the
c subroutine "select_syst" ;
c
c icsort: Fundamental or Observational Catalogue
c 1 = Observational catalogue
c It is assumed that the positions hold for the mean
c epoch of observation as given in the catalogue. Each
c star has thus its individual mean epoch, which may be
c different in right ascension and declination.
c 2 = Fundamental catalogue. In that case it is assumed
c that the positions and porper motions hold for mean
c epoch at the mean equinox of the catalogue, e.g.
c postions and proper motions are given for 1950. and
c they were transferred to 1950 with the proper motions
c as given in the catalogue and with the precession as
c given in the catalogue.
c
c eqk1 : Equinox-Correction in right ascension
c eqk2 : Equinox-Correction for the poper motions in right ascension
c = "motion of the equinox"
c epeqk1: Epoch (and equinox) for which the correction holds.
c ip : Precession, used in the catalogue
c 1 = IAU (1976) System of Astronomical constatns
c 2 = Newcomb's precession
c equcat: Equinox of the catalogue.
c
c
implicit real*8 (a-h,o-z)
c
character*1 answer
c
tra = .261799387799149D+00
trd = .174532925199433D-01
trma = .727220521664304D-04
trmd = .484813681109536D-05
c
c Definition of the catalogue sort:
c =================================
c
write(*,'(//,''Definition of the catalogue sort:'',//,
f ''icsort = 1 : The data are given for the mean epoch of '',
f ''observation and the catalogue equinox'',/,
f '' No proper motions are applied to the '',
f ''positions '',/,
f '' This convention holds for many'',
f '' "observational catalogues" '',//,
f ''icsort = 2 : The data are given for the epoch of the '',
f ''catalogue equinox'',/,
f '' The reduction was made with '',
f ''proper motions as given in the catalogue'',/,
f '' This convention holds e.g. for a '',
f ''fundamental catalogue '',//)')
write(*,'(/,''Please select '')')
read(*,'(i1)') icsort
c
c
c Definition of the various catalogue systems, which can be handled
c =====================================================================
c
c
c Reduction from FK5 system to Hipparcos system
c These parameters are the same for the following 3 subsets of the FK5:
c Basic, Bright Extension + Sup, Faint Extension
c =====================================================================
c
if (isyst.le. 3) then
c
c icsort = 2
ip = 1
equcat = 2000.d0
eqk1 = 0.000d0
eqk2 = 0.000d0
epeqk1 = 2000.d0
c
end if
c
c
c
c Reduction from FK4 system to Hipparcos system
c =============================================
c
if (isyst.eq. 4) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.035d0
eqk2 = 0.085d0
epeqk1 = 1950.d0
c
end if
c
c
c Reduction from FK3 system to Hipparcos system
c =============================================
c
if (isyst.eq. 5) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.035d0
eqk2 = 0.085d0
epeqk1 = 1950.d0
c
end if
c
c
c
c Reduction from NFK system to Hipparcos system
c =============================================
c
if (isyst.eq. 6) then
c
c icsort = 2
ip = 2
equcat = 1875.d0
eqk1 = -0.070d0
eqk2 = +0.085d0
epeqk1 = 1875.d0
c
end if
c
c
c
c Reduction from N30 system to Hipparcos system
c =============================================
c
if (isyst.eq. 7) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.010d0
eqk2 = 0.070d0
epeqk1 = 1930.d0
c
end if
c
c
c
c Reduction from GC system to Hipparcos system
c =============================================
c
if (isyst.eq.8) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.035d0
eqk2 = 0.085d0
epeqk1 = 1950.d0
c
end if
c
c
c
c Reduction from PGC system to Hipparcos system
c ==============================================
c
if (isyst.eq.9) then
c
c icsort = 2
ip = 2
equcat = 1900.d0
eqk1 = -0.050d0
eqk2 = 0.070d0
epeqk1 = 1900.d0
c
end if
c
c
c Reduction from PPM system to Hipparcos system
c ==============================================
c
if (isyst.eq.10) then
c
c icsort = 2
ip = 1
equcat = 2000.d0
eqk1 = 0.000d0
eqk2 = 0.000d0
epeqk1 = 2000.d0
c
end if
c
c
c Reduction from Eichelberger's system to Hipparcos system
c =========================================================
c
if (isyst.eq.11) then
c
c icsort = 2
ip = 2
equcat = 1925.d0
eqk1 = 0.020d0
eqk2 = 0.140d0
epeqk1 = 1925.d0
c
end if
c
c
c Reduction from SAO system to Hipparcos system
c ==============================================
c
if (isyst.eq.12) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.035d0
eqk2 = 0.085d0
epeqk1 = 1950.d0
c
end if
c
c
c Reduction from AGK3 system to Hipparcos system
c ===============================================
c
if (isyst.eq.13) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.035d0
eqk2 = 0.085d0
epeqk1 = 1950.d0
c
end if
c
c
c Reduction from AGK3R system to Hipparcos system
c ================================================
c
if (isyst.eq.14) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
eqk1 = 0.035d0
eqk2 = 0.085d0
epeqk1 = 1950.d0
c
end if
c
c
c Reduction from Perth 70 system (Total catalogue) to Hipparcos system
c =====================================================================
c
if (isyst.eq.15) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
c
eqk1 = +0.050d0
eqk2 = +0.000d0
epeqk1 = 1970.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c Reduction from system of the IRS (on FK4) system to Hipparcos system
c =====================================================================
c
if (isyst.eq.16) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
c
eqk1 = +0.035d0
eqk2 = +0.085d0
epeqk1 = 1950.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c
c Reduction from system of the IRS (on FK5) system to Hipparcos system
c =====================================================================
c
if (isyst.eq.17) then
c
c icsort = 2
ip = 1
equcat = 2000.d0
c
eqk1 = +0.000d0
eqk2 = +0.000d0
epeqk1 = 2000.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c Reduction from system of the system of the AGK2 (printed version)
c to the Hipparcos system
c ==================================================================
c
if (isyst.eq.18) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
c
eqk1 = +0.035d0
eqk2 = +0.085d0
epeqk1 = 1950.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c
c Reduction from system of the system of the AGK2 (revised version)
c to the Hipparcos system
c ==================================================================
c
if (isyst.eq.19) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
c
eqk1 = +0.035d0
eqk2 = +0.085d0
epeqk1 = 1950.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c
c
c Reduction from system of the system of the AGK2A system to the Hipparcos system
c ===============================================================================
c
if (isyst.eq.20) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
c
eqk1 = +0.035d0
eqk2 = +0.085d0
epeqk1 = 1950.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c
c Reduction from system of the ACRS (on FK5) system to Hipparcos system
c =====================================================================
c
if (isyst.eq.21) then
c
c icsort = 2
ip = 1
equcat = 2000.d0
c
eqk1 = +0.000d0
eqk2 = +0.000d0
epeqk1 = 2000.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c Reduction from the SRS system (FK4-based version) to the Hipparcos system
c ==========================================================================
c
if (isyst.eq.22) then
c
c icsort = 2
ip = 2
equcat = 1950.d0
c
eqk1 = +0.035d0
eqk2 = +0.085d0
epeqk1 = 1950.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c
c Reduction from the SRS system (FK5-based version) to the Hipparcos system
c ==========================================================================
c
if (isyst.eq.23) then
c
c icsort = 2
ip = 1
equcat = 2000.d0
c
c eqk1 = +0.035d0
c eqk2 = +0.085d0
epeqk1 = 1950.d0
c
eqk1 = 0.000d0
eqk2 = 0.000d0
c
end if
c
c
c Reduction from the CPC2 system (FK5-based version) to the Hipparcos system
c ==========================================================================
c
if (isyst.eq.24) then
c
c icsort = 2
ip = 1
equcat = 2000.d0
c
eqk1 = +0.000d0
eqk2 = +0.000d0
epeqk1 = 1968.d0
c
c eqk1 = 0.000d0
c eqk2 = 0.000d0
c
end if
c
c
c Read the global catalogue parameters
c
c100 continue
c
write(*,'(''Definition of global catalogue parameters '',//)')
c
c
write(*,'(//,''Defined conventions for the catalogue to '',
f ''be transformed:'',/ )')
c
if (icsort.eq.1) write(*,'(''The catalogue is assumed to be '',
f ''an observational catalogue'',/,
f ''and the positions and proper motion hold for the mean '',
f ''epoch of the individual star.'',/)')
c
if (icsort.eq.2) write(*,'(''The catalogue is assumed to be '',
f ''a fundamental catalogue'',/,
f ''with positions and proper motion given for the epoch '',
f ''of the catalogue equinox,'',/,
f ''and positions reduced to the epoch of the equinox with '',
f ''the proper motions as given in the catalogue.'',/)')
c
if (ip.eq.1) write(*,'(''The precession used in the input '',
f ''catalogue is the IAU (1976) precession'',/,
f ''In this case, it is assumed, that the time unit is the '',
f ''Julian century'',/,
f ''that proper motions provided are per Julian century'',/,
f ''epochs are Julian epochs '',//)')
c
if (ip.eq.2) write(*,'(''The precession used in the '',
f ''input catalogue is Newcomb"s precession.'',/,
f ''In this case, it is assumed, that the time unit is the '',
f ''tropical century'',/,
f ''that proper motions are per tropical century and '',/,
f ''thar epochs are Besselian epochs '',//)')
c
write(*,'(''The catalogue equinox is '',f8.1)') equcat
write(*,'(''Equinox correction in R.A. '',f8.3,'' time sec'',/,
f ''which holds for the epoch '',f8.1)') eqk1,epeqk1
write(*,'(''Motion of the equinox '',f8.3,'' time sec'',
f '' per century '',f8.3)') eqk2
c
write(*,'(//,''Are these definitions correct? y/n'')')
read(*,'(a)') answer
if (answer.eq.'y') goto 100
write(*,'(//,''Make necessary corrections in the '',
f ''subroutine "conventions" '',/,
f ''translate and run the program again '')')
goto 900
100 continue
c
c transform Equinox-Corrections in radian
c
eqk1 = eqk1*trma
eqk2 = eqk2*trma
c
return
c
900 stop
c
end
c
c
subroutine syst_trans (xa,xd,xma,xmd,equ,epa,epd,
f xpar,xvel,ip,sysreg,epverg,
f eqk1,epeqk1,eqk2,icsort,
f ya,yd,yma,ymd)
c
c
c
c Compute and apply the systematic effects for the transition
c from the pre-IAU (1976) system to the IAU (1976) system of constants.
c Zweck: alte Kataloge zunaechst mal global auf das neue
c
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
c
c Unterprogramm hat wohl nicht die hoechste Genauigkeit,
c da ja kleine Restsystemtiken durch den sich anschliessenden
c Katalogvergleich weggebuegelt werden.
c
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
c
c
c Daher wird in der Routine davon ausgegangen, dass die Daten
c des zu transformierenden Katalogs noch den "alten"
c Konventionen entsprechen.
c
c The following effects are taken into account:
c
c 1. Change to IAU (1976) precession
c 2. Elimination of e-terms
c 3. Apply equinox correction (and motion of the equinox)
c 4. Transition to Equinox J2000
c 5. Refer proper motions to on Julian century
c
c Input parameters:
c
c xa, xd: Right ascension and declination of the star in radians
c xma, xmd: Proper motions in radians per century
c epa, epd: Besselian epoch for which xa, xd are given in the catalogue
c equ : Equinox of the catalogue
c epverg : Epoch for which the catalogue comparison holds
c sysreg(1): regional correction in RA at the epoch of the catalogue comparison
c (2): regional correction in Dec at the epoch of the catalogue comparison
c (3): regional correction in RA proper motion
c (4): regional correction in Dec proper motion
c xpar : Parallax in radians
c xvel : Radial velocity in km/s.
c
c eqk1 : Equinox Correction in radian at the epoch epeqk1
c epeqk1 : Besselian epoch for which the Equinox correctin eqk1 holds
c eqk2 : Correction for the "Motion of the Equinox" in rad per century
c ip : Flag describing the Precession used in the catalogue:
c ip = 1 : IAU (1976)
c ip = 2 : Newcomb
c ip = 3 : Struwe
c
c icsort : Kind of catalogue:
c 1 = Observational catalogue. In this case it is assumed that
c the positions hold for the given mean epoch of observation
c of the individual star and no proper motions have been applied.
c 2 = Fundamental catalog. In this case it is assumed that the positions
c have been transformed to the epoch of equinox with the proper motions
c given in the catalogue.
c
c
implicit real*8 (a-h,o-z)
c
real*8 sysreg(4),sysreg2(4)
c
character*12 a12,d12
c
tra = .261799387799149D+00
trd = .174532925199433D-01
trma = .727220521664304D-04
trmd = .484813681109536D-05
trjul = 1.000021356D0
xj2000 = 2451545.D0
c
if (icsort.lt.3) goto 50
write(*,'(//,''ICSORT: Kind of catalogue is wrongly defined,'',/,
f ''The catalogue cannot be reduced'',//)')
goto 900
50 continue
c
c
c In the case of an observational catalogue the positions have to be given
c at the mean epoch of observations. No proper motions are involved.
c In order to make that sure, they are defined here ex[plicitely as zero.
c
if (icsort.ne.1) goto 60
xma = 0.d0
xmd = 0.d0
60 continue
c
c
c Reduce systematic correction to rad and rad/cy
c ==============================================
c
c
cosd = dcos(xd)
sysreg2(1) = sysreg(1)/1000.d0/15.d0/3600.d0*tra/cosd
sysreg2(2) = sysreg(2)/1000.d0/3600.d0*trd
sysreg2(3) = sysreg(3)*100.d0/1000.d0/15.d0/3600.d0*tra/cosd
sysreg2(4) = sysreg(4)*100.d0/1000.d0/3600.d0*trd
c
c
c Elimination of E-Terms
c ======================
c
if (ip.ne.1) then
c
call kstetr (equ,c0,d0,eps)
call eterm (xa,xd,c0,d0,eps,corra,corrd)
xa = xa - corra
xd = xd - corrd
c
end if
c
c
c Treatment of Precession
c =======================
c
c Transformation of the catalogue to Equinox of the epoch of observation
c in the case of an observational catalogue with the precessional values
c as used in the catalogue.
c
c Transformation of epochs in Julina date.
c
if (ip.ne.1) call bejd (epa,epaj)
if (ip.eq.1) call jejd (epa,epaj)
if (ip.ne.1) call bejd (epd,epdj)
if (ip.eq.1) call jejd (epd,epdj)
if (ip.ne.1) call bejd (equ,equj)
if (ip.eq.1) call jejd (equ,equj)
c
c
xp0 = xpar
xv0 = xvel
c
c call radcha(xa,'HH MM SS.SSS',a12,jer,6)
c call radcha(xd,'VDD MM SS.SS',d12,jer,6)
c
call tspv2(ip,equj,epaj,epaj,epaj,xa,xd ,xma, xmd,xp0,xv0,
f xa1,xd1,xma1,xmd1,yp0,yv0,ier)
call tspv2(ip,equj,epdj,epdj,epdj,xa ,xd,xma ,xmd,xp0,xv0,
f xa2,xd2,xma2,xmd2,yp0,yv0,ier)
c
c
if (icsort.ne.1) goto 110
xma1 = 0.d0
xmd1 = 0.d0
xma2 = 0.d0
xmd2 = 0.d0
110 continue
c
c
xp0 = yp0
xv0 = yv0
c
c
c Treatment of Equinox-Correction and global rotation
c ===================================================
c
c Computation of the corrections at the epoch, where the
c catalogue data are effectively given (e.g. 1950 in the case
c of the FK4 catalogue)
c
xda1 = eqk1 + eqk2*(epa-epeqk1)/100.d0
xdma1 = eqk2
xda2 = eqk1 + eqk2*(epd-epeqk1)/100.d0
xdma2 = eqk2
xa1 = xa1 + xda1
xma1 = xma1 + xdma1
xa2 = xa2 + xda2
xma2 = xma2 + xdma2
c
xa1 = xa1 + sysreg2(1) + sysreg2(3)*(epa-epverg)/100.d0
xa2 = xa2 + sysreg2(1) + sysreg2(3)*(epa-epverg)/100.d0
xd1 = xd1 + sysreg2(2) + sysreg2(4)*(epd-epverg)/100.d0
xd2 = xd2 + sysreg2(2) + sysreg2(4)*(epd-epverg)/100.d0
xma1 = xma1 + sysreg2(3)
xma2 = xma2 + sysreg2(3)
xmd1 = xmd1 + sysreg2(4)
xmd2 = xmd2 + sysreg2(4)
c
c
if (icsort.ne.1) goto 120
xma1 = 0.d0
xmd1 = 0.d0
xma2 = 0.d0
xmd2 = 0.d0
120 continue
c
c
xdm1 = 0.d0
xdn1 = 0.d0
xdm2 = 0.d0
xdn2 = 0.d0
c
c
if (icsort.eq.1) goto 300
c
c
c Refer the porper motions to the IAU(1976) precessional rates.
c =============================================================
c
c Case ip = 2 (Newcomb Precession)
c
if (ip.ne.2) goto 200
c
c Unterschied in Praezessionsgeschwindigkeiten definieren
c Vgl. hierzu Unterlagen in der Akte APFS (Programmentwicklung)
c
taua = (epaj-xj2000)/36525.d0
taud = (epdj-xj2000)/36525.d0
c
xdm1 = 1.0368d0 - 0.00176d0*taua - 0.000398d0*taua*taua
xdn1 = 0.4363d0 + 0.00075d0*taua + 0.000153d0*taua*taua
xdm2 = 1.0368d0 - 0.00176d0*taud - 0.000398d0*taud*taud
xdn2 = 0.4363d0 + 0.00075d0*taud + 0.000153d0*taud*taud
c
c
c
c Case ip = 3 (Struwe Precession)
c
200 continue
c
if (ip.ne.3) goto 300
c
c Unterschied in Praezessionsgeschwindigkeiten definieren
c Vgl. hierzu Unterlagen in der Akte ARIGFH (Praezessionsgroessen)
c
taua = (epaj-xj2000)/36525.d0
taud = (epdj-xj2000)/36525.d0
c
xdm1 = 0.4097d0 - 0.05594d0*taua
xdn1 = -0.0759d0 + 0.00972d0*taua
xdm2 = 0.4097d0 - 0.05594d0*taud
xdn2 = -0.0759d0 + 0.00972d0*taud
c
300 continue
c
xdma1 = xdm1 + xdn1*dsin(xa1)*dtan(xd1)
xdmd1 = xdn1*dcos(xa1)
xdma2 = xdm2 + xdn2*dsin(xa2)*dtan(xd2)
xdmd2 = xdn2*dcos(xa2)
c
xma1 = xma1 - xdma1*trmd
xmd1 = xmd1 - xdmd1*trmd
xma2 = xma2 - xdma2*trmd
xmd2 = xmd2 - xdmd2*trmd
c
c
c Refer the porper motions to the Julina century
c ===============================================
c
if (ip.ne.1) then
c
xma1 = xma1*trjul
xmd1 = xmd1*trjul
xma2 = xma2*trjul
xmd2 = xmd2*trjul
c
end if
c
c
if (icsort.ne.1) goto 320
xma1 = 0.d0
xmd1 = 0.d0
xma2 = 0.d0
xmd2 = 0.d0
320 continue
c
c
c Transformation to Equinox J2000 (and epoch J2000 in the case of
c a fundamental catalogue.
c ===============================================================
c
c
jp = 1
c
call tspv2(jp,epaj,epaj,xj2000,xj2000,xa1,xd1,xma1,xmd1,
f xp0,xv0,
f ya1,yd1,yma1,ymd1,
f yp0,yv0,ier)
call tspv2(jp,epdj,epdj,xj2000,xj2000,xa2,xd2,xma2,xmd2,
f xp0,xv0,
f ya2,yd2,yma2,ymd2,
f yp0,yv0,ier)
c
c
if (icsort.ne.1) goto 330
yma1 = 0.d0
ymd2 = 0.d0
330 continue
c
c
ya = ya1
yd = yd2
yma = yma1
ymd = ymd2
c
return
c
c
900 stop
c
end
c
c
DOUBLE PRECISION FUNCTION NORMF(N,M)
C
C BERECHNUNG DER NORMIERUNGSFAKTOREN DER PRODUKTE AUS DEN
C LEGENDRE-POLYNOMEN UND DEN FOURIER-FUNKTIONEN
C
INTEGER N,M
NORMF = DSQRT(2.D0*N + 1.D0)
IF (M.GT.0) NORMF = NORMF*DSQRT(2.D0)
RETURN
END
DOUBLE PRECISION FUNCTION FOUR8 (M,L,ALFA)
C
C FOURIER-FUNKTION (M,L) AN DER STELLE ALFA (IM BOGENMASS) BERECHN.
C
REAL*8 ALFA
INTEGER M,L
FOUR8 = 1.D0
IF (M.EQ.0) GO TO 10
IF (L.EQ.1) FOUR8 = DCOS(L*M*ALFA)
IF (L.EQ.-1) FOUR8 = DSIN(-L*M*ALFA)
10 RETURN
END
REAL FUNCTION XLEG8 (N,X)
C
C BERECHNUNG DES N-TEN LEGENDRE-POLYNOMS AN DER STELLE X
C ACM - ALGORITHMUS NR. 14
C
IMPLICIT REAL*8 (A-Z)
INTEGER I,N,NM1
A = 1.D0
B = X
C = A
IF (N.EQ.1) C = B
IF (N.LT.2) GO TO 20
NM1 = N-1
DO 10 I = 1,NM1
XI = I*1.D0
C = B*X + (XI/(XI+1.D0))*(X*B-A)
A = B
10 B = C
20 XLEG8 = C
RETURN
END
REAL FUNCTION HE*8 (N,X)
C BERECHNUNG DES HERMITE-POLYNOMS H(N) DURCH REKURSION
C
C
REAL*8 A,B,C,X
A = 1.D0
B = X
IF (.NOT. (N.LT.2)) GO TO 10
IF (N.EQ.0) C = A
IF (N.EQ.1) C = B
GO TO 30
10 NM1 = N-1
DO 20 I = 1,NM1
C = X*B - I*A
A = B
20 B = C
30 HE = C
RETURN
END
REAL FUNCTION NORMH*8 (N)
C BERECHNUNG VON 1/N - FAKULTAET
REAL*8 F
F = 1.D0
IF (N.LT.1) GO TO 20
DO 10 I = 1,N
10 F = F*DSQRT(I*1.D0)
20 NORMH = 1.D0/F
RETURN
END
c
SUBROUTINE INDUX(K,P,N,M,L)
C AUS K WERDEN P,N,M,L BERECHNET
C
INTEGER P
IN = K
L = -1
IF (IN.GT.0) L = 1
IN = IN*L
P = IN/1 000 000
N = (IN-P*1 000 000)/1 000
M = (IN-P*1 000 000 - N*1 000)
RETURN
END
c
c
SUBROUTINE ROTINF(ALPHA,DELTA,XOM1,XOM2,XOM3,XDA,XDD)
C
C BERECHNUNG DER KOORDINATENAENDERUNGEN BEI EINER INFINI-
C TESIMALEN ROTATION DER SPHAERE UM DIE INFINITESIMALEN WINKEL
C XOM1, XOM2, XOM3.
C ORIENTIERUNG WIE BEI LINDEGREN (NDAC/LO/148, 5. JULI 1993)
C
C
C GILT ALSO NUR FUER KLEINE WINKEL XOM1, XOM2, XOM3 !!!!!!!!!!!
C =============
C
C EINGABE: ALPHA, DELTA = SPHAERISCHE KOORDIANTEN (DES STERNS)
C IN RADIAN.
C XOM1 = A23 = INFINITESIMALE ROTATION UM X1-ACHSE
C XOM2 = A13 = INFINITESIMALE ROTATION UM X2-ACHSE
C XOM3 = A12 = INFINITESIMALE ROTATION UM X3-ACHSE
C AUSGABE: XDA, XDD = INFINITESIMALE KOORDIANTENAENDERUNG AN DER
C STELLE ALPHA, DELTA IN DEN GLEICHEN EINHEITEN
C WIE DIE ROTATIONSWINKEL
C
IMPLICIT REAL*8 (A-H,O-Z)
C
XSA = DSIN(ALPHA)
XCA = DCOS(ALPHA)
XSD = DSIN(DELTA)
XCD = DCOS(DELTA)
XDA = (-XOM3*XCD + XOM2*XSA*XSD + XOM1*XCA*XSD)/XCD
XDD = ( + XOM2*XCA - XOM1*XSA)
RETURN
END
c
C *******************************************************************
C ***** *
C ***** UNTERPROGRAMM CHARAD *
C ***** *
C ***** UMWANDLUNG EINES WINKELS, DER IN GRAD ODER STUNDEN ODER *
C ***** MINUTEN (ZEIT ODER BOGEN-) ODER IN SEKUNDEN (ZEIT ODER *
C ***** BOGEN-) ODER IN EINER SINNVOLLEN KOMBINATION DER GENANN- *
C ***** TEN EINHEITEN ALS CHARCTERSTRING ANGEGEBEN IST, IN *
C ***** BOGENMASS (DOPPELTGENAUE REALGROESSE). *
C ***** *
C *******************************************************************
C ***** *
C ***** AUFGERUFEN WIRD DAS UNTERPROGRAM MIT *
C ***** *
C ***** CALL CHARAD (C1,C2,RDG,IER,LUNERR) *
C ***** *
C ***** WOBEI BEDEUTEN : *
C ***** C1 = WINKEL IN D,H,M,S ETC C*(*) EINGABE *
C ***** C2 = FORMAT VON C1 C*(*) EINGABE *
C ***** RDG = WINKEL IM BOGENMASS R*8 AUSGABE *
C ***** IER = FEHLERMELDUNG I3 AUSGABE *
C ***** LUNERR = LOGICAL UNIT NUMBER FUER DIE AUSGABE EINER *
C ***** FEHLERMELDUNG I2 EINGABE *
C ***** *
C *******************************************************************
SUBROUTINE CHARAD(CHA,FORMAT,RAD,IER,LUNERR)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DATA FZ,SZ /15.D0,60.D0/
CHARACTER*1 CH(20),KCH(20)
CHARACTER*(*) CHA,FORMAT,ierstr*3
CHARACTER*70 text0,text1,text2,text3,text4,text5,text6,text7
CHARACTER*70 text8
ILE1= LEN(CHA)
ILE2= LEN(FORMAT)
C
text0 = 'von charad behoben, Fehlerangabe nur zur Information'
text1 = 'XSTAR8(3 bzw 6) auf -9999.d0 KENNX8(3 bzw 6) auf 0 setz
*en'
text2 = 'korrigierter Wert ist in die Datei korrektur einzutragen'
text3 = 'Koordinate fehlt'
text4 = 'XSTAR8(3 bzw 6) auf -8888.d0 KENNX8(3 bzw 6) auf 0 setz
*en'
text5 = 'falls es zuviele sind charad ueberspringen'
text6 = 'Fehler in der Formatangabe des Charad-Aufrufes'
text7 = 'im Leseprogramm zu korrigieren'
text8 = 'ergaenzte Koordinate ist in die Datei korrektur einzutrag
*en'
C Erstmal den Formatstring auf Grossbuchstaben bringen
do i=1,ile2
call upcase(format(i:i))
end do
IF(ILE1.NE.ILE2) THEN
IER = 93
GOTO 900
END IF
DO 99 I=ILE1,1,-1
IF (CHA(I:I).NE.' ') GOTO 9
99 CONTINUE
IER = 99
GOTO 900
9 LCHA=I
C *******************************************************************
C ***** *
C ***** RA = 1 WENN DER WINKEL IN GRAD, BOGENMINUTEN ODER *
C ***** -SEKUNDEN ANGEBEN IST *
C ***** = 15 WENN DER WINKEL IN STUNDEN, ZEITMINUTEN *
C ***** -SEKUNDEN ANGEBEN IST *
C ***** = 0 WENN EIN FEHLER AUFGETRETEN IST, DER MIT *
C ***** IER-MELDUNG ENDET *
C ***** ZU BEGINN DES PROGRAMMS WIRD RA=1 GESETZT *
C ***** *
C ***** IER = 0 KEIN FEHLER BEI DER EINGABE *
C ***** > 0 FEHLER BEI DER EINGABE *
C ***** WENN MOD(10) = 0 FEHLER IN FORMATEINAGE (C2) *
C ***** WENN MOD(10) > 0 FEHLER IN WERTEINGABE (C1) *
C ***** ZU BEGINN DES PROGRAMMS WIRD IER=0 GESETZT *
C ***** *
C ***** IVA = 1 POSITIVES VORZEICHEN (' ' , + , & ) *
C ***** = -1 NEGATIVES VORZEICHEN ( - ) *
C ***** ZU BEGINN DES PROGRAMMS WIRD IVA=1 GESETZT *
C ***** *
C ***** IVK = 0 KEINE UEBERLOCHUNG DER LETZTEN STELLE IN C1 *
C ***** = 1 UEBERLOCHUNG DER LETZTEN STELLE IN C1 *
C ***** ZU BEGINN DES PROGRAMMS WIRD IVK=0 GESET
C ***** *
C *******************************************************************
RA=1.D0
PI = 3.1415926535897932D0
IER = 0
IVA=1
IVK=0
A=3600.D0
RADZ=0.D0
DO 10 I=1,LCHA
CH(I)=CHA(I:I)
KCH(I)=FORMAT(I:I)
10 CONTINUE
C *******************************************************************
C ***** *
C ***** ABFRAGE DES FORMATS, IN WELCHER EINHEIT C1 ANGEGEBEN IST *
C ***** *
C ***** H = STUNDEN (REKTASZENSION) *
C ***** D = GRAD (DEKLINATION, GAL.BREITE ETC.) *
C ***** L = LAENGE (GAL.LAENGE) *
C ***** T = ZEIT- (-MINUTEN ODER -SEKUNDEN) *
C ***** A = BOGEN (-MINUTEN ODER -SEKUNDEN) *
C ***** V = VORZEICHEN (POS: ' ','+','&') *
C ***** (NEG: '-') *
C ***** X = UEBERLOCHUNG (DER 1. ODER LETZTEN STELLE) *
C ***** M = MINUTEN (ZUSAMMEN MIT T ODER HH : ZEIT
C ***** (MIT A, DD ODER ' ' : BOGENMINUTEN) *
C ***** S = SEKUNDEN (ZUSAMMEN MIT T ODER HH : ZEITSEKUNDEN) *
C ***** (MIT A, DD ODER ' ' : BOGENSEKUNDEN) *
C ***** *
C *******************************************************************
IF (KCH(2).EQ.'H') GOTO 100
IF (KCH(2).EQ.'D') GOTO 200
IF (KCH(2).EQ.'L') GOTO 300
J=1
IF (KCH(1).NE.'T') GOTO 20
RA=FZ
GOTO 50
20 IF (KCH(1).EQ.'A') GOTO 50
IF (KCH(1).EQ.'V') GOTO 30
IF (KCH(1).EQ.'X') GOTO 50
IF (KCH(LCHA).EQ.'X') GOTO 40
IF (KCH(1).EQ.'Y') GOTO 61
IF (KCH(LCHA).EQ.'Y') GOTO 60
GOTO 50
61 RA=FZ
GOTO 50
60 RA=FZ
GOTO 40
29 IF (KCH(1).EQ.'D') GOTO 401
IF (KCH(1).EQ.'M') GOTO 411
IF (KCH(1).EQ.'S') GOTO 421
GOTO 90
50 IF (KCH(2).EQ.'M') GOTO 210
A=SZ
IF (KCH(2).EQ.'S') GOTO 220
90 IER = 90
GOTO 900
30 IVA=1
IF (CH(1).EQ.'-') IVA=-1
J=J+1
IF (KCH(2).EQ.'T') RA=FZ
IF (KCH(3).EQ.'M') GOTO 210
IF (KCH(3).EQ.'S') GOTO 220
GOTO 90
40 IVK=1
CALL INPUT(CH(LCHA),Z,IV,IE)
IF (IE.EQ.0) GOTO 41
IF (IE.GT.0) IER = 91
IF (IE.LT.0) IER = 92
GOTO 900
41 B=Z
IVA=IV
IF (LCHA.EQ.2) GOTO 29
GOTO 50
C *******************************************************************
C ***** *
C ***** UMRECHNUNG VON STUNDENWINKELN *
C ***** *
C *******************************************************************
100 RA=FZ
J=0
ZW=0.D0
DO 102 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 1010
IF (IE.GT.0) IER = 151
IF (IE.LT.0) IER = 152
GOTO 900
1010 ZZ = Z * DBLE(K)
ZW = ZW + ZZ
IF (ZW.LE.24.) GOTO 101
IER = 111
GOTO 900
101 RADZ=RADZ + ZZ * 3600.D0
IF (J.GE.LCHA) GOTO 900
102 CONTINUE
103 J=J+1
IF (KCH(J).EQ.' ') GOTO 103
IF (KCH(J).EQ.'M') GOTO 110
IF (KCH(J).EQ.'.') GOTO 150
IF (KCH(J).EQ.'H') GOTO 151
IER = 110
GOTO 900
110 J=J-1
ZW = 0.D0
DO 112 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 1100
IF (IE.GT.0) IER = 151
IF (IE.LT.0) IER = 152
GOTO 900
1100 ZZ = Z * DBLE(K)
ZW = ZW + ZZ
IF (ZW.LE.SZ) GOTO 111
IER = 121
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
111 RADZ = RADZ + ZZ * SZ
IF (J.GE.LCHA) GOTO 900
112 CONTINUE
A=A/SZ
115 J=J+1
IF (KCH(J).EQ.' ') GOTO 115
IF (KCH(J).EQ.'S') GOTO 120
IF (KCH(J).EQ.'.') GOTO 150
IF (KCH(J).EQ.'M') GOTO 151
IER = 120
GOTO 900
120 J=J-1
ZW = 0.D0
DO 122 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 1200
IF (IE.GT.0) IER = 151
IF (IE.LT.0) IER = 152
GOTO 900
1200 ZZ = Z * DBLE (K)
ZW = ZW + ZZ
IF (ZW.LE.SZ) GOTO 121
IER = 131
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
121 RADZ = RADZ + Z * DBLE(K)
IF (J.GE.LCHA) GOTO 900
122 CONTINUE
A=A/SZ
123 J=J+1
IF (KCH(J).EQ.'.') GOTO 150
IF (KCH(J).EQ.'S') GOTO 151
IF (KCH(J).EQ.' ') GOTO 123
IER = 130
GOTO 900
150 J=J+1
151 CONTINUE
DO 152 K=J,LCHA
CALL INPUT (CH(K),Z,IV,IE)
IF (IE.EQ.0) GOTO 1520
IF (IE.GT.0) IER = 91
IF (IE.LT.0) IER = 92
GOTO 900
1520 RADZ = RADZ +A * Z/(10.D0**DBLE(K-J+1))
152 CONTINUE
GOTO 900
C *******************************************************************
C ***** *
C ***** UMRECHNUNG VON GRAD *
C ***** *
C *******************************************************************
200 IVK=0
IV =0
J=0
IVA=1
N=0
IF (KCH(LCHA).EQ.'X') GOTO 270
IF (KCH(1).EQ.'X') GOTO 203
IF (KCH(1).EQ.'D') GOTO 203
IF (CH(1).EQ.'-') GOTO 201
IF (CH(1).EQ.'&') GOTO 202
IF (CH(1).EQ.'+') GOTO 202
IF (CH(1).EQ.' ') GOTO 202
IER = 200
GOTO 900
201 IVA=-1
202 J=J+1
203 IF (KCH(2).EQ.'X') GOTO 401
ZW = 0.D0
DO 204 K=10,1,-9
J=J+1
N=N+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 2203
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
2203 IF(IVA.LT.0) GOTO 2204
IF(K.EQ.10) IVA=IV
2204 ZZ = Z * DBLE(K)
ZW = ZW + ZZ
IF (ZW.LE.90.D0) GOTO 204
IER = 211
GOTO 900
204 RADZ = (RADZ + ZZ * 3600.D0)
IF (J.GE.LCHA) GOTO 900
205 CONTINUE
206 J=J+1
IF (KCH(J).EQ.' ') GOTO 206
IF (KCH(J).EQ.'M') GOTO 210
IF (KCH(J).EQ.'.') GOTO 250
IF (KCH(J).EQ.'D') GOTO 251
IF (KCH(J).EQ.'X') GOTO 251
IER = 210
GOTO 900
210 IF (KCH(J+1).EQ.'X') GOTO 411
J=J-1
ZW = 0.D0
DO 212 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 2110
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
2110 ZZ = Z * DBLE(K)
ZW = ZW + ZZ
IF (ZW.LE.SZ) GOTO 211
IER = 221
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
211 RADZ = RADZ + Z * DBLE(K) * SZ
IF (J.GE.LCHA) GOTO 900
212 CONTINUE
A=A/SZ
215 J=J+1
IF (KCH(J).EQ.' ') GOTO 215
IF (KCH(J).EQ.'S') GOTO 220
IF (KCH(J).EQ.'.') GOTO 250
IF (KCH(J).EQ.'M') GOTO 251
IF (KCH(J).EQ.'X') GOTO 251
IER = 220
GOTO 900
220 IF (KCH(J+1).EQ.'X') GOTO 421
J=J-1
ZW = 0.D0
DO 222 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 2200
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
2200 ZZ = Z * DBLE(K)
ZW = ZW + ZZ
IF (ZW.LE.SZ) GOTO 221
IER = 231
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
221 RADZ = RADZ + Z * DBLE(K)
IF (J.GE.LCHA) GOTO 900
222 CONTINUE
A=A/SZ
223 J=J+1
IF(KCH(J).EQ.'.') GOTO 250
IF(KCH(J).EQ.'S') GOTO 251
IF(KCH(J).EQ.' ') GOTO 223
IF(KCH(J).EQ.'X') GOTO 251
IER = 230
GOTO 900
250 J=J+1
251 J1 = LCHA
IF (IVK.NE.1) GOTO 252
J1 = J1-1
RADZ = RADZ + A * B/(10.D0**DBLE(LCHA-J+1))
252 IF (J.GT.J1) GOTO 900
DO 253 K=J,J1
CALL INPUT (CH(K),Z,IV,IE)
IF (IE.EQ.0) GOTO 2520
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
2520 RADZ = RADZ + A * Z/(10.D0** DBLE(K-J+1))
253 CONTINUE
GOTO 900
270 IVK=1
CALL INPUT (CH(LCHA),Z,IV,IE)
IF (IE.EQ.0) GOTO 2700
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
2700 B=Z
IVA=IV
GOTO 203
401 CALL INPUT(CH(1),Z,IV,IE)
IF (IE.EQ.0) GOTO 4010
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
4010 ZW=10.D0*Z
IF (ZW.LE.90.D0) GOTO 402
IER = 211
GOTO 900
402 RADZ = RADZ + Z * 36000.D0 + B * 3600.D0
GOTO 900
411 CALL INPUT(CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 4110
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
4110 ZW=10.D0*Z
IF (ZW.LE.SZ) GOTO 412
IER = 221
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
412 RADZ = RADZ + Z * 600.D0 + B * SZ
GOTO 900
421 CALL INPUT(CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 4210
IF (IE.GT.0) IER = 251
IF (IE.LT.0) IER = 252
GOTO 900
4210 ZW=10.D0*Z
IF (ZW.LE.SZ) GOTO 422
IER = 231
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C go to 900
422 RADZ = RADZ + Z * 10.D0 + B
GOTO 900
C *******************************************************************
C ***** *
C ***** UMRECHNUNG VON GAL.LAENGE *
C ***** *
C *******************************************************************
300 J=1
IVA=1
IF(CH(1).EQ.'-') IVA=-1
ZW=0.D0
DO 302 K=2,0,-1
J=J+1
CALL INPUT(CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 3000
IF (IE.GT.0) IER = 351
IF (IE.LT.0) IER = 352
GOTO 900
3000 ZZ=Z*10.D0**DBLE(K)
ZW=ZW+ZZ
IF (ZZ.LE.360.D0) GOTO 301
IER = 311
GOTO 900
301 RADZ=RADZ+ZZ*3600.D0
IF (J.GE.LCHA) GOTO 900
302 CONTINUE
303 J=J+1
IF (KCH(J).EQ.' ') GOTO 303
IF (KCH(J).EQ.'M') GOTO 310
IF (KCH(J).EQ.'.') GOTO 350
IF (KCH(J).EQ.'L') GOTO 351
IER = 310
GOTO 900
310 J=J-1
ZW = 0.D0
DO 312 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 3100
IF (IE.GT.0) IER = 351
IF (IE.LT.0) IER = 352
GOTO 900
3100 ZZ = Z * DBLE(K)
ZW = ZW + ZZ
IF (ZW.LE.SZ) GOTO 311
IER = 321
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
311 RADZ = RADZ + ZZ * SZ
IF (J.GE.LCHA) GOTO 900
312 CONTINUE
A=A/SZ
315 J=J+1
IF (KCH(J).EQ.' ') GOTO 315
IF (KCH(J).EQ.'S') GOTO 320
IF (KCH(J).EQ.'.') GOTO 350
IF (KCH(J).EQ.'M') GOTO 351
IER = 320
GOTO 900
320 J=J-1
ZW = 0.D0
DO 322 K=10,1,-9
J=J+1
CALL INPUT (CH(J),Z,IV,IE)
IF (IE.EQ.0) GOTO 3200
IF (IE.GT.0) IER = 351
IF (IE.LT.0) IER = 352
GOTO 900
3200 ZZ = Z * DBLE (K)
ZW = ZW + ZZ
IF (ZW.LE.SZ) GOTO 321
IER = 331
WRITE (LUNERR,1111) CHA,IER,FORMAT
C
print *, text0
C
C GOTO 900
321 RADZ = RADZ + Z * DBLE(K)
IF (J.GE.LCHA) GOTO 900
322 CONTINUE
A=A/SZ
323 J=J+1
IF (KCH(J).EQ.'.') GOTO 350
IF (KCH(J).EQ.'S') GOTO 351
IF (KCH(J).EQ.' ') GOTO 323
IER = 330
GOTO 900
350 J=J+1
351 CONTINUE
DO 352 K=J,LCHA
CALL INPUT (CH(K),Z,IV,IE)
IF (IE.EQ.0) GOTO 3510
IF (IE.GT.0) IER = 351
IF (IE.LT.0) IER = 352
GOTO 900
3510 RADZ = RADZ +A * Z/(10.D0**DBLE(K-J+1))
352 CONTINUE
GOTO 900
900 CONTINUE
IF (IER.EQ.0) GOTO 910
IF (IER.EQ.121) GOTO 910
IF (IER.EQ.131) GOTO 910
IF (IER.EQ.221) GOTO 910
IF (IER.EQ.231) GOTO 910
WRITE (LUNERR,1111) CHA,IER,FORMAT
1111 FORMAT (//, 'FEHLER AUFGETRETEN BEI UMRECHNUNG VON'/
*1X,A,' IN RAD. FEHLERMELDUNG = ',I3/1X,A)
C
write(ierstr,'(i3)')ier
C
if(ier.eq.99)then
print *, text3
print *, text4
print *, text8
print *, text5
end if
C
if(ierstr(3:3).eq.'0')then
print *, text6
print *, text7
else
if(ier.gt.100)then
print *, text1
print *, text2
end if
end if
C
RAD = 0.D0
GOTO 920
910 RAD = ( RADZ / 3600.D0 ) * RA * ( PI / 180.D0)
IF (IVA.EQ.-1) RAD=-RAD
920 CONTINUE
RETURN
END
C *******************************************************************
C ***** *
C ***** UNTERPROGRAMM : INPUT *
C ***** *
C ***** UMWANDLUNG EINER CHARCTERGROESSE IN EINE REALGROESSE *
C ***** *
C *******************************************************************
C ***** *
C ***** AUFGERUFEN WIRD DAS UNTERPROGRAMM MIT *
C ***** *
C ***** CALL INPUT(C,R,IV,IE) *
C ***** *
C ***** WOBEI BEDEUTEN *
C ***** C = CHARACTERGROESSE A*1 EINGABE *
C ***** R = REALGROESSE (DOPPELTGENAU) R*8 AUSGABE *
C ***** IV = VORZEICHEN A*1 AUSGABE *
C ***** WIRD DURCH INPUT NUR GEAENDERT, WENN EINE *
C ***** ENTSPRECHENDE UEBERLOCHUNG VORLIEGT *
C ***** IE = FEHLERMELDUNG I1 AUSGABE *
C ***** = 1 WENN C NICHT MIT DEM VORGEGEBENEN *
C ***** ZEICHENSATZ UEBEREINSTIMMT *
C ***** BEWIRKT, DASS IER = 301 *
C ***** = -1 WENN C = '.' (PUNKT) *
C ***** BEWIRKT, DASS IER = 311 *
C ***** *
C *******************************************************************
SUBROUTINE INPUT(C,Z,IV,IE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*1 C,IZI(10),IBUP(10),IBUN(10)
DATA IZI/'0','1','2','3','4','5','6','7','8','9'/
DATA IBUN/ '}','J','K','L','M','N','O','P','Q','R'/
DATA IBUP/ '{','A','B','C','D','E','F','G','H','I'/
IV =1
IE = 0
IF (C.NE.' ') GOTO 9
I=1
GOTO 21
9 CONTINUE
DO 10 I=1,10
IF (C.EQ.IZI(I)) GOTO 21
IF (C.EQ.IBUP(I)) GOTO 21
IF (C.EQ.IBUN(I)) GOTO 20
10 CONTINUE
IE = 1
IF (C.EQ.'.') IE=-1
Z = 0.D0
GOTO 22
20 IV=-1
21 Z = DBLE(I-1)
22 CONTINUE
RETURN
END
c
C********************************************************************
SUBROUTINE RADCHA(RADX,C2,CWERT,IER,LUNERR)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*(*) C2,CWERT
CHARACTER*10 CBT
CHARACTER*1 CHR(20),CZ(4)
CHARACTER*2 CH,CM,CS,CD,CV
DATA CZ/' ','1','2','3'/
ILE1= LEN(C2)
ILE2= LEN(CWERT)
C Erstmal den Formatstring auf Grossbuchstaben bringen
do i=1,ile1
call upcase(c2(i:i))
end do
RAD=RADX
IE1=0
INULF=0
DO 99 I=ILE1,1,-1
IF (C2(I:I).NE.' ') GOTO 98
99 CONTINUE
IER=90
GOTO 900
98 LAE=I
IF (C2(I:I).NE.'F') GOTO 97
LAE=I-1
INULF=1
97 CONTINUE
IF (C2(ILE1:ILE1).EQ.'F') ILE1 = ILE1 - 1
IF(ILE1.NE.ILE2) THEN
IER = 93
GOTO 900
END IF
DO 96 I=1,ILE2
CWERT(I:I)=' '
CHR(I)=' '
96 CONTINUE
DO 95 I=1,LAE
CHR(I)=C2(I:I)
95 CONTINUE
RA=1.
PI = 3.1415926535897932D0
IER=0
IVA=1
IH=0
IRUND=0
K1=0
IF (RAD.LT.0.D0) IVA=-1
RAD = RAD * IVA
IF (CHR(2).EQ.'H') IH=1
IF (CHR(1).EQ.'T') IH=1
IF (CHR(2).EQ.'T') IH=1
IF (IH.EQ.1) GOTO 50
C********************************************************************
C
C
C
C********************************************************************
R1=RAD*180.D0/PI
RD=DINT(R1)
RDB=R1-RD
R2=RDB*60.D0
RM=DINT(R2)
RMB=R2-RM
R3=RMB*60.D0
RS=DINT(R3)
RSB=R3-RS
IF (RS.LT.60.D0) GOTO 41
RS=RS-60.D0
RM=RM+1.D0
41 IF (RM.LT.60.D0) GOTO 42
RM=RM-60.D0
RD=RD+1.D0
42 IF (CHR(2).EQ.'M') GOTO 60
IF (CHR(2).EQ.'A') GOTO 60
IF (CHR(2).EQ.'L' .AND. RD.LT.360.D0) GOTO 60
IF (CHR(2).EQ.'D' .AND. RD.LT.90.D0) GOTO 60
IF (CHR(2).EQ.'L') RD=RD-360.D0
IF (CHR(2).EQ.'D') RD=RD-90.D0
GOTO 42
C********************************************************************
C
C
C
C********************************************************************
50 R1=RAD*180.D0/(PI*15.D0)
RH=DINT(R1)
RHB=R1-RH
R2=RHB*60.D0
RM=DINT(R2)
RMB=R2-RM
R3=RMB*60.D0
RS=DINT(R3)
RSB=R3-RS
IF (RS.LT.60.D0) GOTO 51
RS=RS-60.D0
RM=RM+1.D0
51 IF (RM.LT.60.D0) GOTO 52
RM=RM-60.D0
RH=RH+1.D0
52 IF (RH.LT.24.D0) GOTO 60
RH=RH-24.D0
GOTO 52
C********************************************************************
C
C
C
C********************************************************************
60 IF (CHR(1).EQ.'H') GOTO 100
IF (CHR(2).EQ.'D') GOTO 200
IF (CHR(2).EQ.'L') GOTO 500
IF (CHR(2).EQ.'T') GOTO 600
IF (CHR(2).EQ.'A') GOTO 700
IF (CHR(2).EQ.'M') GOTO 700
C********************************************************************
C
C
C
C********************************************************************
100 IE1=0
IF (INULF.EQ.1) GOTO 180
CALL INPUT1(RH,CH,IE)
IE1=IE1+IE
CALL INPUT1(RM,CM,IE)
IE1=IE1+IE
CALL INPUT1(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 190
IER=101
GOTO 900
180 CALL INPUT2(RH,CH,IE)
IE1=IE1+IE
CALL INPUT2(RM,CM,IE)
IE1=IE1+IE
CALL INPUT2(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 190
IER=102
GOTO 900
190 CWERT(1:2)=CH
J=2
101 J=J+1
IF (J.GT.LAE) GOTO 401
IF (CHR(J).NE.' ') GOTO 102
CWERT(J:J)=' '
GOTO 101
102 IF (CHR(J).NE.'M') GOTO 110
C CWERT=CWERT//CM
J1=J
J2=J+1
CWERT(J1:J2)=CM
J=J+1
103 J=J+1
IF (J.GT.LAE) GOTO 402
IF (CHR(J).NE.' ') GOTO 104
CWERT(J:J)=' '
GOTO 103
104 IF (CHR(J).NE.'S') GOTO 120
J1=J
J2=J+1
CWERT(J1:J2)=CS
J=J+1
C J=J+2
IF (J.GT.LAE) GOTO 403
105 J=J+1
IF (J.GT.LAE) GOTO 403
IF (CHR(J).NE.' ') GOTO 130
CWERT(J:J)=' '
GOTO 105
110 IF (CHR(J).NE.'.') GOTO 112
CWERT(J:J)='.'
111 J=J+1
IF (J.GT.LAE) GOTO 401
IF (CHR(J).NE.' ') GOTO 112
CWERT(J:J)=' '
GOTO 111
112 IF (CHR(J).NE.'H') GOTO 113
BT=RHB
K1=1
GOTO 350
113 IF (CHR(J).NE.' ') GOTO 114
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 401
GOTO 112
114 IER=111
GOTO 900
120 IF (CHR(J).NE.'.') GOTO 122
CWERT(J:J)='.'
121 J=J+1
IF (J.GT.LAE) GOTO 402
IF (CHR(J).NE.' ') GOTO 122
CWERT(J:J)=' '
GOTO 121
122 IF (CHR(J).NE.'M') GOTO 123
BT=RMB
K1=2
GOTO 350
123 IF (CHR(J).NE.' ') GOTO 124
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 402
GOTO 122
124 IER=121
GOTO 900
130 IF (CHR(J).NE.'.') GOTO 132
CWERT(J:J)='.'
131 J=J+1
IF (J.GT.LAE) GOTO 403
IF (CHR(J).NE.' ') GOTO 132
CWERT(J:J)=' '
GOTO 131
132 IF (CHR(J).NE.'S') GOTO 133
BT=RSB
K1=3
GOTO 350
133 IF (CHR(J).NE.' ') GOTO 134
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 403
GOTO 132
134 IER=131
GOTO 900
C********************************************************************
C
C
C
C********************************************************************
200 IE1=0
IF (INULF.EQ.1) GOTO 280
CALL INPUT1(RD,CD,IE)
IE1=IE1+IE
CALL INPUT1(RM,CM,IE)
IE1=IE1+IE
CALL INPUT1(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 290
IER=201
GOTO 900
280 CALL INPUT2(RD,CD,IE)
IE1=IE1+IE
CALL INPUT2(RM,CM,IE)
IE1=IE1+IE
CALL INPUT2(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 290
IER=202
GOTO 900
290 CWERT(2:3)=CD
J=3
201 J=J+1
IF (J.GT.LAE) GOTO 411
IF (CHR(J).NE.' ') GOTO 202
CWERT(J:J)=' '
GOTO 201
202 IF (CHR(J).NE.'M') GOTO 210
J1=J
J2=J+1
CWERT(J1:J2)=CM
J=J+1
203 J=J+1
IF (J.GT.LAE) GOTO 412
IF (CHR(J).NE.' ') GOTO 204
CWERT(J:J)=' '
GOTO 203
204 IF (CHR(J).NE.'S') GOTO 220
J1=J
J2=J+1
CWERT(J1:J2)=CS
J=J+1
IF (J.GT.LAE) GOTO 413
205 J=J+1
IF (J.GT.LAE) GOTO 413
IF (CHR(J).NE.' ') GOTO 230
CWERT(J:J)=' '
GOTO 205
210 IF (CHR(J).NE.'.') GOTO 212
CWERT(J:J)='.'
211 J=J+1
IF (J.GT.LAE) GOTO 411
IF (CHR(J).NE.' ') GOTO 212
CWERT(J:J)=' '
GOTO 211
212 IF (CHR(J).NE.'D') GOTO 213
BT=RDB
K1=4
GOTO 350
213 IF (CHR(J).NE.' ') GOTO 214
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 411
GOTO 212
214 IER=211
GOTO 900
220 IF (CHR(J).NE.'.') GOTO 222
CWERT(J:J)='.'
221 J=J+1
IF (J.GT.LAE) GOTO 412
IF (CHR(J).NE.' ') GOTO 222
CWERT(J:J)=' '
GOTO 221
222 IF (CHR(J).NE.'M') GOTO 223
BT=RMB
K1=5
GOTO 350
223 IF (CHR(J).NE.' ') GOTO 224
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 412
GOTO 222
224 IER=221
GOTO 900
230 IF (CHR(J).NE.'.') GOTO 232
CWERT(J:J)='.'
231 J=J+1
IF (J.GT.LAE) GOTO 413
IF (CHR(J).NE.' ') GOTO 232
CWERT(J:J)=' '
GOTO 231
232 IF (CHR(J).NE.'S') GOTO 233
BT=RSB
K1=6
GOTO 350
233 IF (CHR(J).NE.' ') GOTO 234
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 413
GOTO 232
234 IER=231
GOTO 900
C********************************************************************
C
C
C
C********************************************************************
500 IE1=0
IF (RD.LT.360.D0) GOTO 591
RD=RD-360.D0
GOTO 500
591 RD1=DINT(RD/100.D0)
RD2=RD-RD1*100.D0
IF (INULF.EQ.1) GOTO 580
CALL INPUT1(RD2,CD,IE)
IE1=IE1+IE
CALL INPUT1(RM,CM,IE)
IE1=IE1+IE
CALL INPUT1(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 590
IER=301
GOTO 900
580 CALL INPUT2(RD2,CD,IE)
IE1=IE1+IE
CALL INPUT2(RM,CM,IE)
IE1=IE1+IE
CALL INPUT2(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 590
IER=302
GOTO 900
590 CONTINUE
DO 599 I=0,3
IF (RD1.EQ.DBLE(I)) GOTO 598
599 CONTINUE
IER=301
GOTO 900
598 CWERT(2:2)=CZ(I+1)
CWERT(3:4)=CD
IF (CWERT(2:2).NE.' '.AND.CWERT(3:3).EQ.' ') CWERT(3:3)='0'
IF (INULF.NE.1) GOTO 597
IF (CWERT(2:2).EQ.' ') CWERT(2:2)='0'
IF (CWERT(3:3).EQ.' ') CWERT(3:3)='0'
597 J=4
501 J=J+1
IF (J.GT.LAE) GOTO 421
IF (CHR(J).NE.' ') GOTO 502
CWERT(J:J)=' '
GOTO 501
502 IF (CHR(J).NE.'M') GOTO 510
J1=J
J2=J+1
CWERT(J1:J2)=CM
J=J+1
503 J=J+1
IF (J.GT.LAE) GOTO 422
IF (CHR(J).NE.' ') GOTO 504
CWERT(J:J)=' '
GOTO 503
504 IF (CHR(J).NE.'S') GOTO 520
J1=J
J2=J+1
CWERT(J1:J2)=CS
J=J+1
IF (J.GT.LAE) GOTO 423
505 J=J+1
IF (J.GT.LAE) GOTO 423
IF (CHR(J).NE.' ') GOTO 530
CWERT(J:J)=' '
GOTO 505
510 IF (CHR(J).NE.'.') GOTO 512
CWERT(J:J)='.'
511 J=J+1
IF (J.GT.LAE) GOTO 421
IF (CHR(J).NE.' ') GOTO 512
CWERT(J:J)=' '
GOTO 511
512 IF (CHR(J).NE.'L') GOTO 513
BT=RDB
K1=7
GOTO 350
513 IF (CHR(J).NE.' ') GOTO 514
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 421
GOTO 512
514 IER=311
GOTO 900
520 IF (CHR(J).NE.'.') GOTO 522
CWERT(J:J)='.'
521 J=J+1
IF (J.GT.LAE) GOTO 422
IF (CHR(J).NE.' ') GOTO 522
CWERT(J:J)=' '
GOTO 521
522 IF (CHR(J).NE.'M') GOTO 523
BT=RMB
K1=8
GOTO 350
523 IF (CHR(J).NE.' ') GOTO 524
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 422
GOTO 522
524 IER=321
GOTO 900
530 IF (CHR(J).NE.'.') GOTO 532
CWERT(J:J)='.'
531 J=J+1
IF (J.GT.LAE) GOTO 423
IF (CHR(J).NE.' ') GOTO 532
CWERT(J:J)=' '
GOTO 531
532 IF (CHR(J).NE.'S') GOTO 533
BT=RSB
K1=9
GOTO 350
533 IF (CHR(J).NE.' ') GOTO 534
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 423
GOTO 532
534 IER=331
GOTO 900
C********************************************************************
C
C
C
C********************************************************************
600 IE1=0
IF (INULF.EQ.1) GOTO 680
CALL INPUT1(RM,CM,IE)
IE1=IE1+IE
CALL INPUT1(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 690
IER=101
GOTO 900
680 CALL INPUT2(RM,CM,IE)
IE1=IE1+IE
CALL INPUT2(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 690
IER=102
GOTO 900
690 CWERT(2:3)=CM
J=3
603 J=J+1
IF (J.GT.LAE) GOTO 432
IF (CHR(J).NE.' ') GOTO 604
CWERT(J:J)=' '
GOTO 603
604 IF (CHR(J).NE.'S') GOTO 620
J1=J
J2=J+1
CWERT(J1:J2)=CS
J=J+1
IF (J.GT.LAE) GOTO 433
605 J=J+1
IF (J.GT.LAE) GOTO 433
IF (CHR(J).NE.' ') GOTO 630
CWERT(J:J)=' '
GOTO 605
620 IF (CHR(J).NE.'.') GOTO 622
CWERT(J:J)='.'
621 J=J+1
IF (J.GT.LAE) GOTO 432
IF (CHR(J).NE.' ') GOTO 622
CWERT(J:J)=' '
GOTO 621
622 IF (CHR(J).NE.'M') GOTO 623
BT=RMB
K1=11
GOTO 350
623 IF (CHR(J).NE.' ') GOTO 624
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 432
GOTO 622
624 IER=121
GOTO 900
630 IF (CHR(J).NE.'.') GOTO 632
CWERT(J:J)='.'
631 J=J+1
IF (J.GT.LAE) GOTO 433
IF (CHR(J).NE.' ') GOTO 632
CWERT(J:J)=' '
GOTO 631
632 IF (CHR(J).NE.'S') GOTO 633
BT=RSB
K1=12
GOTO 350
633 IF (CHR(J).NE.' ') GOTO 634
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 433
GOTO 632
634 IER=131
GOTO 900
C********************************************************************
C
C
C
C
C********************************************************************
700 IE1=0
IF (INULF.EQ.1) GOTO 780
CALL INPUT1(RM,CM,IE)
IE1=IE1+IE
CALL INPUT1(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 790
IER=201
GOTO 900
780 CALL INPUT2(RM,CM,IE)
IE1=IE1+IE
CALL INPUT2(RS,CS,IE)
IE1=IE1+IE
IF (IE1.EQ.0) GOTO 790
IER=202
GOTO 900
790 CWERT(2:3)=CM
J=3
703 J=J+1
IF (J.GT.LAE) GOTO 442
IF (CHR(J).NE.' ') GOTO 704
CWERT(J:J)=' '
GOTO 703
704 IF (CHR(J).NE.'S') GOTO 720
J1=J
J2=J+1
CWERT(J1:J2)=CS
J=J+1
IF (J.GT.LAE) GOTO 443
705 J=J+1
IF (J.GT.LAE) GOTO 443
IF (CHR(J).NE.' ') GOTO 730
CWERT(J:J)=' '
GOTO 705
720 IF (CHR(J).NE.'.') GOTO 722
CWERT(J:J)='.'
721 J=J+1
IF (J.GT.LAE) GOTO 442
IF (CHR(J).NE.' ') GOTO 722
CWERT(J:J)=' '
GOTO 721
722 IF (CHR(J).NE.'M') GOTO 723
BT=RMB
K1=13
GOTO 350
723 IF (CHR(J).NE.' ') GOTO 724
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 442
GOTO 722
724 IER=221
GOTO 900
730 IF (CHR(J).NE.'.') GOTO 732
CWERT(J:J)='.'
731 J=J+1
IF (J.GT.LAE) GOTO 443
IF (CHR(J).NE.' ') GOTO 732
CWERT(J:J)=' '
GOTO 731
732 IF (CHR(J).NE.'S') GOTO 733
BT=RSB
K1=15
GOTO 350
733 IF (CHR(J).NE.' ') GOTO 734
CWERT(J:J)=' '
J=J+1
IF (J.GT.LAE) GOTO 443
GOTO 732
734 IER=231
GOTO 900
C********************************************************************
C
C
C
C********************************************************************
350 IF (IRUND.NE.0) GOTO 380
K=LAE-J+1
A=DBLE(K)
BT=DINT(BT*10.**A+0.5D0)
AT=BT/(10.D0**A)
IF (AT.LT.1.D0) GOTO 380
BT=BT-10.D0**(A+1.D0)
IRUND=1
IF (K1.EQ.1) GOTO 361
IF (K1.EQ.4) GOTO 364
IF (K1.EQ.7) GOTO 367
IF (K1.EQ.2) GOTO 362
IF (K1.EQ.5) GOTO 365
IF (K1.EQ.8) GOTO 368
IF (K1.EQ.3) GOTO 363
IF (K1.EQ.6) GOTO 366
IF (K1.EQ.9) GOTO 369
IF (K1.EQ.11) GOTO 371
IF (K1.EQ.12) GOTO 372
IF (K1.EQ.14) GOTO 374
IF (K1.EQ.15) GOTO 375
361 RH=RH+1
IF (RH.EQ.24.) RH=0.D0
IF (RH.EQ.25.) RH=1.D0
GOTO 100
362 RM=RM+1
IF (RM.LT.60.D0) GOTO 100
RM=RM-60.D0
GOTO 361
363 RS=RS+1
IF (RS.LT.60.D0) GOTO 100
RS=RS-60.D0
GOTO 362
364 RD=RD+1
GOTO 200
365 RM=RM+1
IF (RM.LT.60.D0) GOTO 200
RM=RM-60.D0
GOTO 364
366 RS=RS+1
IF (RS.LT.60.D0) GOTO 200
RS=RS-60.D0
GOTO 365
367 RD=RD+1
IF (RD.GE.360.D0) RD=RD-360.D0
GOTO 500
368 RM=RM+1
IF (RM.LT.60.D0) GOTO 500
RM=RM-60.D0
GOTO 367
369 RS=RS+1
IF (RS.LT.60.D0) GOTO 500
RS=RS-60.D0
GOTO 368
371 RM=RM+1
IF (RM.LT.60.D0) GOTO 600
RM=RM-60.D0
GOTO 600
372 RS=RS+1
IF (RS.LT.60.D0) GOTO 600
RS=RS-60.D0
GOTO 371
374 RM=RM+1
IF (RM.LT.60.D0) GOTO 700
RM=RM-60.D0
GOTO 700
375 RS=RS+1
IF (RS.LT.60.D0) GOTO 700
RS=RS-60.D0
GOTO 374
C********************************************************************
C
C
C
C********************************************************************
380 CALL INPUT3(BT,K,CBT,IE)
J1=J
J2=LAE
CWERT(J1:J2)=CBT
GOTO 900
C********************************************************************
C
C
C
C********************************************************************
401 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RHB1=RHB+0.5D0
IF (RHB1.LT.1.0D0) GOTO 900
GOTO 361
402 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RMB1=RMB+0.5D0
IF (RMB1.LT.1.0D0) GOTO 900
GOTO 362
403 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RSB1=RSB+0.5D0
IF (RSB1.LT.1.0D0) GOTO 900
GOTO 363
411 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RDB1=RDB+0.5D0
IF (RDB1.LT.1.0D0) GOTO 900
GOTO 364
412 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RMB1=RMB+0.5D0
IF (RMB1.LT.1.0D0) GOTO 900
GOTO 365
413 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RSB1=RSB+0.5D0
IF (RSB1.LT.1.0D0) GOTO 900
GOTO 366
421 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RDB1=RDB+0.5D0
IF (RDB1.LT.1.0D0) GOTO 900
GOTO 367
422 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RMB1=RMB+0.5D0
IF (RMB1.LT.1.0D0) GOTO 900
GOTO 368
423 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RSB1=RSB+0.5D0
IF (RSB1.LT.1.0D0) GOTO 900
GOTO 369
432 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RMB1=RMB+0.5D0
IF (RMB1.LT.1.0D0) GOTO 900
GOTO 371
433 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RSB1=RSB+0.5D0
IF (RSB1.LT.1.0D0) GOTO 900
GOTO 372
442 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RMB1=RMB+0.5D0
IF (RMB1.LT.1.0D0) GOTO 900
GOTO 374
443 IF (IRUND.EQ.1) GOTO 900
IRUND=1
RSB1=RSB+0.5D0
IF (RSB1.LT.1.0D0) GOTO 900
GOTO 375
C********************************************************************
C
C
C
C********************************************************************
900 CONTINUE
IF (IER.EQ.0) GOTO 910
WRITE(LUNERR,1111) RAD,C2,IER
1111 FORMAT (1X, 'FEHLER AUFGETRETEN BEI UMRECHNUNG VON'/
*1X,F11.8,' IN ',A/1X,'FEHLERMELDUNG = ',I3//)
RAD = 0.D0
GOTO 920
910 CONTINUE
CV='+'
IF (CHR(1).EQ.'N') CV=' '
IF (IVA.EQ.-1) CV='-'
IF (IH.NE.1) CWERT(1:1)=CV
IF (IH.EQ.1 .AND. CHR(2).EQ.'T') CWERT(1:1)=CV
920 CONTINUE
999 CONTINUE
RETURN
END
C********************************************************************
C
C
C
C********************************************************************
SUBROUTINE INPUT1(Z,C,IE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*2 IZI(100),C
DATA IZI/' 0',' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9',
*'10','11','12','13','14','15','16','17','18','19',
*'20','21','22','23','24','25','26','27','28','29',
*'30','31','32','33','34','35','36','37','38','39',
*'40','41','42','43','44','45','46','47','48','49',
*'50','51','52','53','54','55','56','57','58','59',
*'60','61','62','63','64','65','66','67','68','69',
*'70','71','72','73','74','75','76','77','78','79',
*'80','81','82','83','84','85','86','87','88','89',
*'90','91','92','93','94','95','96','97','98','99'/
IE = 0
DO 10 I=0,99
IF (Z.EQ.DBLE(I)) C=IZI(I+1)
10 CONTINUE
RETURN
END
C********************************************************************
C
C
C
C********************************************************************
SUBROUTINE INPUT2(Z,C,IE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*2 IZI(100),C
DATA IZI/'00','01','02','03','04','05','06','07','08','09',
*'10','11','12','13','14','15','16','17','18','19',
*'20','21','22','23','24','25','26','27','28','29',
*'30','31','32','33','34','35','36','37','38','39',
*'40','41','42','43','44','45','46','47','48','49',
*'50','51','52','53','54','55','56','57','58','59',
*'60','61','62','63','64','65','66','67','68','69',
*'70','71','72','73','74','75','76','77','78','79',
*'80','81','82','83','84','85','86','87','88','89',
*'90','91','92','93','94','95','96','97','98','99'/
IE = 0
DO 10 I=0,99
IF (Z.EQ.DBLE(I)) C=IZI(I+1)
10 CONTINUE
RETURN
END
C********************************************************************
C
C
C********************************************************************
SUBROUTINE INPUT3(B,L,CC,IE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION Z(9)
CHARACTER*10 CC
CHARACTER*1 IZ(10),IZI(10)
DATA IZ/'0','1','2','3','4','5','6','7','8','9'/
DO 99 I=1,10
CC(I:I)=' '
99 CONTINUE
IF (L.GT.9) GOTO 90
DO 10 I=1,9
Z(I)=0
IZI(I)=' '
10 CONTINUE
GOTO(1,2,3,4,5,6,7,8,9) L
9 Z(9)=DINT(B/10.**8.D0)
8 Z(8)=DINT((B-Z(9)*10.D0**8.D0)/10.**7.D0)
7 Z(7)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0)/10.**6.D0)
6 Z(6)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0-
F Z(7)*10.D0**6.D0)/10.**5.D0)
5 Z(5)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0-
F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0)/10.**4.D0)
4 Z(4)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0-
F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0-
F Z(5)*10.D0**4.D0)/1000.)
3 Z(3)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0-
F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0-
F Z(5)*10.D0**4.D0-Z(4)*10.D0**3.D0)/100.)
2 Z(2)=DINT((B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0-
F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0-
F Z(5)*10.D0**4.D0-Z(4)*10.D0**3.D0-
F Z(3)*10.D0**2.D0)/10.)
1 Z(1)= DINT(B-Z(9)*10.D0**8.D0-Z(8)*10.D0**7.D0-
F Z(7)*10.D0**6.D0-Z(6)*10.D0**5.D0-
F Z(5)*10.D0**4.D0-Z(4)*10.D0**3.D0-
F Z(3)*10.D0**2.D0-Z(2)*10.D0**1.D0)
DO 12 I=0,9
DO 11 J=1,9
IF (Z(J).EQ.DBLE(I)) IZI(J)=IZ(I+1)
11 CONTINUE
12 CONTINUE
CC(1:1)=IZI(L)
DO 13 I=L-1,1,-1
L1=L-I+1
CC(L1:L1)=IZI(I)
13 CONTINUE
GOTO 14
90 IE=1
14 CONTINUE
RETURN
C********************************************************************
C
C
C
C********************************************************************
END
C Wandeln von Klein- nach Grossbuchstaben
Subroutine Upcase(c)
character c*1
integer i
i = ichar(c)
if ((i.ge.97).and.(i.le.122)) then
c = char(i-32)
end if
end
c
SUBROUTINE JEJD(JE,JD)
C
C UMWANDLUNG DER JULIANISCHEN EPOCHE IN JULIANISCHES DATUM
C
REAL*8 JE,JD
C
JD = (JE-2000.D0)*365.25D0 + 2451545.D0
C
RETURN
END
c
c
C 1. Name: Subroutine BEJD U. Bastian, Okt. 1984
C
C 2. Zweck: Besselsche Epochen in Julianische Daten umwandeln.
C
C 3. Aufruf:
C CALL BEJD(BE,TJD)
C
C 4. Parameter:
C
C BE Input Besselsche Epoche (Jahre) Real*8
C
C TJD Output Julianisches (Ephemeris) Datum (Tage) Real*8
C
C
SUBROUTINE BEJD(BE,TJD)
IMPLICIT REAL*8 (A-Z)
C
TJD=(BE-1900.D0)*365.242198781D0+2415020.31352D0
RETURN
END
c
SUBROUTINE TSPV2 (IP,XEQ,XEP,YEQ,YEP,XA,XD,XM,XMS,XP0,XVEL0,
1 YA,YD,YM,YMS,YP0,YVEL0,IER)
C
C EPOCHEN- UND AEQUINOX-UEBERTRAGUNG VON OERTERN U. EIGENBEWEGUNGEN
C
C EOCHEN-UEBERTRAGUNG MIT BENUTZUNG DER RAUMGESCHWINDIGKEIT
C..............................................................
C ...................
C
C IP : STEUERT, WELCHE PRAEZESSIONSWERTE GENOMMEN WERDEN.
C IP= 1 : NEUE IAU(1976) WERTE;
C IP = 2 : NEWCOMBS PRAEZESSION;
C IP = 3 : STRUVE'S PRAEZESSION;
C
C EINGABE: XEQ = JULIANISCHES DATUM DES AUSGANGS-AEQUINOXES
C YEQ = JULIANISCHES DATUM DES END-AEQUINOXES
C XEP = JULIANISCHES DATUM DER AUSGANGS-EPOCHE
C YEP = JULIANISCHES DATUM DER END-EPOCHE
C XA = REKTASZENSION FOER EPOCHE XEP, AEQUINOX XEQ
C XD = DEKLINATION FUER EPOCHE XEP, AEQUINOX XEQ
C XM = E.B. IN ALPHA FUER EPOCHE YEP, AEQUINOX XEQ
C XMS = E.B. IN DELTA FUER EPOCHE XEP, AEQUINOX XEQ
C XP0 = PARALLAXE ZUR EPOCHE XEP IM BOGENMASS
C XVEL0 = RADIALGESCHWINDIGKEIT ZUR EPOCH XEP
C
C AUSGABE: YA = REKTASZENSION FUER EPOCHE YEP, AEQUINOX YEQ
C YD = DEKLINATION FUER EPOCHE YEP, AEQUINOX YEQ
C YM = E.B. IN ALPHA FUER EPOCHE YEP, AEQUINOX YEQ
C YMS = E.B. IN DELTA FUER EPOCHE YEP, AEQUINOX YEQ
C YP0 = PARALLAXE ZUR EPOCHE XEP IM BOGENMASS
C YVEL0 = RADIALGESCHWINDIGKEIT ZUR EPOCH XEP
C RETURN-CODE IER: = 0: RAD.VEL. UND PAR. UNGLEICH O
C = 1: RAD.VEL. ODER PAR. GLEICH 0
C
C ALLE WINKEL IM BOGENMASS
C ZEITEINHEIT = JULIANISCHES JAHRHUNDERT
C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
REAL*8 PX(3,3)
C
C
C
C EPOCHENUEBERTRAGUNG DER OERTER UND EIGENBEWEGUNGEN
C
CALL TPMSPC (XEP,YEP,XA ,XD ,XM ,XMS ,XP0,XVEL0,
1 XA2,XD2,XM2,XMS2,YP0,YVEL0,IER)
C
C PRAEZESSIONSUEBERTRAGUNG VON OERTERN UND EIGENBEWEGUNGEN
C
CALL KPRZ (IP,XEQ,YEQ,ZET,TETA,ZETA)
CALL PRZMTX( ZET,TETA,ZETA,PX)
C
CALL PRECES(XA2,XD2,XM2,XMS2,PX,XA3,XD3,XM3,XMS3)
C
YA = XA3
YD = XD3
YM = XM3
YMS=XMS3
RETURN
END
c
SUBROUTINE KSTETR(BY,C0,D0,EPS)
C BESTIMMUNG DES "ELLIPTISCHEN TEILS DER ABERRATION DAY NUMBERS"
C EPS WIRD BEZOGEN AUF DAS MITTLERE AEQUNOX FUER BEGINN DES
C BESSELSCHEN JAHRES DER EPOCHE TA. DIE ZEIT IST IN
C EPHEMERIDENZEIT (JULIANISCH) ZU RECHNEN (EXPL. SUPPL. P.98, 489)
C BY = BESSELSCHES JAHR UND TEILE DAVON, FUER DAS GERECHNET WIRD
C DT = DIFFERENZ ZWISCHEN 1900.0 UND BY IN JULIANISCHEN JAHRHUNDERTEN
C DTEPS = DIFFERENZ ZWISCHEN 1900.0 UND BY IN TROPISCHEN JAHRHUNDERTEN
C KAPPA = ABERRATIONS-KONSTANTE = 20.496" (AB 1968)
C
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 T,BY,DT,DTEPS,PERL,EXC,EPS,C0,D0,KAPPA
T = 3.1415926535897932D0/(180.D0*3600.D0)
KAPPA = 20.496D0
DT = 365.2422D0*(BY-1900.D0)/36525.D0
DTEPS = (BY - 1900.D0)/100.D0
PERL = 281.22083 3D0 + 1.71917 5D0*DT + 0.00045 2778D0*DT*DT +
1 0.00000 33333D0*DT*DT*DT
EXC = 0.01675 104D0 - 0.00004 180D0*DT - 0.00000 0126D0*DT*DT
EPS = 23.45229 4D0 - 0.01301 25D0*DTEPS -
1 0.00000 164D0*DTEPS*DTEPS + 0.00000 0503D0*DTEPS*DTEPS*DTEPS
PERL = PERL - 180.D0
PERL = PERL*3600.D0*T
EPS = EPS*3600.D0*T
KAPPA = KAPPA*T
C0 = KAPPA*EXC*DCOS(PERL)*DCOS(EPS)
D0 = KAPPA*EXC*DSIN(PERL)
RETURN
END
c
SUBROUTINE ETERM(ALF,DEL,C0,D0,EPS,C1,C2)
C
C BERECHNUNG DER KORREKTIONEN, DIE VON DER ELLIPTIZITAET
C DER ERDBAHN HERRUEHREN (LINEARER TEIL)
C
C LITERATUR: WOOLARD & CLEMENCE, SPHER. ASTR., P. 114
C EXPL. SUPPL. P. 48,151
C EPS = SCHIEFE DER EKLIPTIK, WIRD IN GESONDERTEM U.P. BERECHNET
C C0,D0 = "ELLIPTISCHER TEIL DER ABERRATION DAY NUMBERS", WERDEN
C KSTETR BERECHNET
C C,CS,D,DS = STAR NUMBERS
C A,D = ORT IM BOGENMASS
C C1,C2 = E-TERM DER ELLIPTISCHEN ABERRATION (LINEARER TEIL)
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 ALF,DEL,C1,C2,C0,D0,C,CS,D,DS,EPS
C = DCOS(ALF)/DCOS(DEL)
D = DSIN(ALF)/DCOS(DEL)
CS = DTAN(EPS)*DCOS(DEL) - DSIN(ALF)*DSIN(DEL)
DS = DCOS(ALF)*DSIN(DEL)
C1 = C*C0 + D*D0
C2 = CS*C0 + DS*D0
RETURN
END
c
SUBROUTINE TPMSPC(T0,T1,A0,D0,XMY0,XMYS0,PAR0,VEL0,
1 A1,D1,XMY1,XMYS1,PAR1,VEL1,IER)
C
C
C EPOCHENUEBERTRAGUNG VON OERTERN UND EIGENBEWEGUNGEN UNTER DER
C ANNAHME KONSTANTER RAUMGESCHWINDIGKEIT.
C
C EINGABEDATEN:
C..................
C T0, T1: ANFANGS-UND ENDEPOCHE IN JULIANISCHEM DATUM
C A0, D0: ORT ZUR ANFANGSEPOCHE, IM BOGENMASS
C XMY0,XMYS0: EIG.BEW.KOMP. IM BOGENMASS, PRO JUL. JAHRH.
C PAR: PARALLAXE IM BOGENMASS ZUR ANFANGS-EPOCHE
C VEL0: RADIALGSCHW. IN KM/SEC ZUR ANFANGSEPOCHE
C
C AUSGABEDATEN:
C..................
C WIE EINGABEDATEN, MIT '1' STATT '0'
C RETURN-CODE IER: = 0: RAD.VEL. UND PAR. UNGLEICH 0
C.....................
C = 1: RAD.VEL. ODER PAR. GLEICH 0
C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION R0(3),R0P(3),XA(3),XD(3),R(3)
C
PI = 3.141592653589793D0
TR = PI/180.D0
TR1 = 3600.D0/TR
TR2 = 1.D0/TR1
IER = 0
F1 = 21.094953D0
F2 = 4.8481368D-6
1700 FORMAT(/,1X,3D25.10)
C
C WRITE(6,1700) T0,T1
C WRITE(6,1700) A0,D0
C WRITE(6,1700) XMY0,XMYS0
X1 = DABS(PAR0)
X2 = DABS(VEL0)
IF ( ( X1 .LT.1.D-15) .OR. ( X2 .LT.1.D-15) ) IER = 1
C
C
C EIGENBEWEGUNGEN UND PARALLAXE IN BOGENSEKUNDEN UMRECHNEN.
C
AA = A0
DA = D0
XMYA = XMY0
XMYSA = XMYS0
PARA = PAR0
VELA = VEL0
XMYA = XMYA*TR1
XMYSA=XMYSA*TR1
PARA = PARA*TR1
C WRITE(6,1700) XMYA,XMYSA,PARA
C
C
C VEKTOR ZUM STERN UND GESCHWINDIGKEITSVEKTOR BERECHNEN,
C REDUZIERT AUF DIE EINHEITSSPHAERE.
C
CALL DIRCOS(AA,DA,R0)
C WRITE(6,1720) AA,DA
C WRITE(6,1720) R0
1720 FORMAT(/,1X,3D25.13)
CALL DLDA(AA,DA,XA)
CALL DLDD(AA,DA,XD)
C WRITE(6,1720) XA
C WRITE(6,1720) XD
DO 20 I = 1,3
20 R0P(I) = XMYA*XA(I)+XMYSA*XD(I) + F1*PARA*VELA*R0(I)
C WRITE(6,1720) R0P
C
C
C EPOCHENDIFFERENZ IN JULIANISCHEN JAHRHUNDERTEN:
C
DT = (T1-T0)/36525.D0
C WRITE(6,1720) DT
C
C
C VEKTOR ZUM STERN ZUR ENDEPOCHE
C
DO 80 I = 1,3
80 R(I) = R0(I) + R0P(I)*DT*F2
C WRITE(6,1720) R
C
C
C ENTFERNUNG DES STERNS IM BOGENMASS ZUR ENDEPOCHE:
C
CALL SCLPRD(R,R,E,3)
E = DSQRT(E)
C WRITE(6,1720) E
PAR1 = PARA/E
C WRITE(6,1720) PAR1
C
C
C UEBERGANG ZUM EINHEITSVEKTOR
C ALPHA , DELTA AUR ENDEPOCHE BERECHNEN.
C
DO 100 I = 1,3
100 R(I) = R(I)/E
C WRITE(6,1700) R
C
CALL ANGLE(R,A1,D1)
C WRITE(6,1700) A1,D1
C
C
C NEUE RADIALGESCHWINDIGKEIT ZUR ENDEPOCHE:
C
CALL SCLPRD(R0P,R,VEL1,3)
IF (IER.EQ.0) VEL1 = VEL1/(F1*PAR1)
IF (IER.EQ.1) VEL1 = VELA
C WRITE(6,1700) VEL1
C
C
C NEUE EIGENBEW.KOMP. ZUR ENDEPOCHE BERECHNEN:
C
CD = DCOS(D1)
CALL DLDA(A1,D1,XA)
CALL DLDD(A1,D1,XD)
CALL SCLPRD(R0P,XA,XMY1,3)
CALL SCLPRD(R0P,XD,XMYS1,3)
C WRITE(6,1700) XMY1,XMYS1
XMY1 = XMY1/(E*CD*CD)
XMYS1 = XMYS1/E
C WRITE(6,1700) XMY1,XMYS1
C
C EIGENBEW. UND PARALLAXE IN BOGENMASS UMWANDELN.
C
XMY1 = XMY1*TR2
XMYS1 = XMYS1*TR2
PAR1 = PAR1*TR2
C
RETURN
END
c
SUBROUTINE KPRZ(IP,T0,T1,ZET,TETA,ZETA)
C
C BESTIMMUNG DER PRAEZESSIONSWINKEL NACH VERSCHIEDENEN KONSTANTEN
C JE NACH WAHL DES PARAMETERS IP:
C
C IP = 1 : NEUE IAU(1976) PRAEZESSION,
C IP = 2 : NEWCOMB'S PRAEZESSION,
C IP = 3 : STRUVE'S PRAEZESSION.
C
C EINGABE: T0 = JULIANISCHES DATUM DER AUSGANGSEPOCHE
C T1 = JULIANISCHES DATUM DER ENDEPOCHE
C AUSGABE: PRAEZESSIONSWINKEL IM BOGENMASS
C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
PI = 3.141592653589793D0
C
C EPOCHENDIFFERENZEN IN JULIANISCHE JAHRHUNDERTE UMRECHNEN
C
XJ2000 = 2451545.D0
DT0 = (T0-XJ2000)/36525.D0
DT1 = (T1-T0)/36525.D0
C
C................................................................
IF (IP.NE.1) GOTO 200
C
C BESTIMMUNG DER PRAEZESSINSWINKEL NACH DEN ANGABEN
C BEI LIESKE, ASTRON.&ASTROPHYS., 73, 282-284,(1979).
C ......................................................
C
ZETA = ( 2306.2181D0 + 1.39656D0*DT0 - 0.000139D0*DT0*DT0)*DT1 +
1 ( 0.30188D0 - 0.000344D0*DT0)*DT1*DT1 +
2 0.017998D0*DT1*DT1*DT1
ZET = ( 2306.2181D0 + 1.39656D0*DT0 - 0.000139D0*DT0*DT0)*DT1 +
1 ( 1.09468D0 + 0.000066D0*DT0)*DT1*DT1 +
2 0.018203D0*DT1*DT1*DT1
TETA = ( 2004.3109D0 - 0.85330D0*DT0 - 0.000217D0*DT0*DT0)*DT1 +
1 (-0.42665D0 - 0.000217D0*DT0)*DT1*DT1 -
2 0.041833D0*DT1*DT1*DT1
C
GOTO 700
C
200 CONTINUE
C
C ....................................................
IF (IP.NE.2) GOTO 300
C
C PRAEZESSIONSWINKEL NACH NEWCOMBS PRAEZESSIONSWERTEN
C ....................................................
C
ZETA = ( 2305.69970D0 + 1.3974397D0*DT0 + 0.000060D0*DT0*DT0)*DT1+
1 ( 0.3020079D0 - 0.000270D0*DT0)*DT1*DT1 +
2 0.0179962D0*DT1*DT1*DT1
ZET = ( 2305.69970D0 + 1.3974397D0*DT0 + 0.000060D0*DT0*DT0)*DT1+
1 ( 1.095432D0 + 0.00039D0*DT0)*DT1*DT1 +
2 0.018326D0*DT1*DT1*DT1
TETA = ( 2003.874600D0 - 0.8540465D0*DT0 - 0.000370D0*DT0*DT0)*DT1
1 +(-0.427073D0 - 0.000370D0*DT0)*DT1*DT1 -
2 0.041803D0*DT1*DT1*DT1
C
GOTO 700
C
300 CONTINUE
C
C ...........................................................
IF (IP.NE.3) GOTO 400
C
C BESTIMMUNG DER PRAEZESSIONSWINKEL NACH STRUVE'S PRAEZESSION
C ...........................................................
C
ZETA = ( 2306.013272D0 + 1.424560853 D0*DT0)*DT1 +
1 0.316313512D0*DT1*DT1
ZET = ( 2306.013272D0 + 1.424560853 D0*DT0)*DT1 +
1 1.108347345D0*DT1*DT1
TETA = ( 2004.386800D0 - 0.8630368665D0*DT0)*DT1 +
1 (-0.4315184332D0)*DT1*DT1
C
400 CONTINUE
C
700 CONTINUE
C
C UMWANDLUNG INS BOGENMASS
C
X = (180.D0*3600.D0)
ZETA = (ZETA*PI)/X
TETA = (TETA*PI)/X
ZET = (ZET *PI)/X
RETURN
END
c
SUBROUTINE PRZMTX (ZET,TETA,ZETA,PREC)
C
C BERECHNUNG DER PRAEZESSIONSMATRIX
C
C EINGABE DER PRAEZESSINSWINKEL IM BOGENMASS
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 PREC(3,3)
C
SZET = DSIN(ZET)
CZET = DCOS(ZET)
SZETA = DSIN(ZETA)
CZETA = DCOS(ZETA)
STETA = DSIN(TETA)
CTETA = DCOS(TETA)
C
PREC(1,1) = +CZETA*CTETA*CZET - SZETA*SZET
PREC(1,2) = -SZETA*CTETA*CZET - CZETA*SZET
PREC(1,3) = -STETA*CZET
PREC(2,1) = +CZETA*CTETA*SZET + SZETA*CZET
PREC(2,2) = -SZETA*CTETA*SZET + CZETA*CZET
PREC(2,3) = -STETA*SZET
PREC(3,1) = +CZETA*STETA
PREC(3,2) = -SZETA*STETA
PREC(3,3) = +CTETA
RETURN
END
c
SUBROUTINE PRECES (XA,XD,XMUE,XMUES,PX,YA,YD,YMUE,YMUES)
C
C PRAEZESSIONSUEBERTRAGUNG VON EIGENBEWEGUNGEN UND OERTERN
C EINE PRAEZESSION DER EIGENBEWEGUNGEN ERFORDERT STETS AUCH DIE
C PRAEZESSION DER OERTER.
C
C EINGABE: XA = ALPHA FUER AUSGANGSAEQUINOX
C XD = DELTA FUER AUSGANGSAEQUINOX
C XMUE = E.B. IN ALPHA FUER AUSGANGSAEQUINOX
C XMUES = E.B. IN DELTA FUER AUSGANGSAEQUINOX
C PX = PRAEZESSIONSMATRIX
C AUSGABE: YA = ALPHA FUER ENDAEQUINOX
C YD = DELTA FUER ENDAEQUINOX
C YMUE = E.B. IN ALPHA FUER ENDAEQUINOX
C YMUES = E.B. IN DELTA FUER ENDAEQUINOX
C ALLE WINKEL IM BOGENMASS
C
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 PX(3,3),VA(3,1),VD(3,1),XRES(3,1),YRES(3,1),ZRES(3,1)
C
C PRAEZESSION DES ORTES
C
CALL PRCPOS (XA,XD,YA,YD,PX)
C
C BERECHNUNG DER VEKTOREN DER ABGELEITETEN RICHTUNGSKOSINUSSE
C
CALL DLDA (XA,XD,VA)
CALL DLDD (XA,XD,VD)
C
C PRAEZESSION DER EIGENBEWEGUNGEN
C
CALL GMPRD (PX,VA,XRES,3,3,1)
CALL SMPY (XRES,XMUE,YRES,3,1)
CALL GMPRD (PX,VD,XRES,3,3,1)
CALL SMPY (XRES,XMUES,ZRES,3,1)
CALL GMADD (YRES,ZRES,XRES,3,1)
C
C AUFLOESUNG NACH DEN EIGENBEWEGUNGS-KOMPONENETEN
C
CALL DLDA (YA,YD,VA)
CALL SCLPRD (XRES,VA,YMUE,3)
C
YMUE = YMUE/(DCOS(YD)*DCOS(YD))
C
CALL DLDD (YA,YD,VD)
CALL SCLPRD (XRES,VD,YMUES,3)
C
YMUES = YMUES
RETURN
END
c
C
SUBROUTINE DIRCOS (A,D,V)
C
C BERECHNUNG DES VEKTORS V DER RICHTUNGSKOSINUSSE
C
C EINGABE: A, D = ALPHA, DELTA IM BOGENMASS
C AUSGABE: V = VEKTOR DER RICHTUNGSKOSINUSSE
C
IMPLICIT REAL*8 (A-H,O-Z)
C
REAL*8 V(3,1)
C
C
V(1,1) = DCOS(A)*DCOS(D)
V(2,1) = DSIN(A)*DCOS(D)
V(3,1) = DSIN(D)
RETURN
END
c
C
SUBROUTINE DLDA (A,D,V)
C
C BERECHNUNG DES NACH ALPHA ABGELEITETEN VEKTORS DER RICHTUNGSKOS.
C
C EINGABE : A, D = ALPHA, DELTA (IM BOGENMASS)
C AUSGABE : V = VEKTOR DER ABGELETETEN RICHTUNGSKOSINUSSE
C
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 V(3,1)
C
V(1,1) = -DCOS(D)*DSIN(A)
V(2,1) = +DCOS(D)*DCOS(A)
V(3,1) = 0.D0
RETURN
END
c
C
SUBROUTINE DLDD (A,D,V)
C
C BERECHNUNG DES NACH DELTA ABGELEITETEN VEKTORS DER
C RICHTUNGSKOSINUSSE
C
C EINGABE: A, D = ALPHA,DELTA (IM BOGENMASS)
C AUSGABE: V = VEKTOR DER RICHTUNGSKOSINUSSE
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 V(3,1)
C
V(1,1) = -DSIN(D)*DCOS(A)
V(2,1) = -DSIN(D)*DSIN(A)
V(3,1) = +DCOS(D)
RETURN
END
c
C
SUBROUTINE SCLPRD (A,B,R,M)
C
C SKALRAPRODUKT R = A.B ZWEIER VEKTOREN
C
C A = NAME DES 1. EINGABEVEKTORS (M-MAL-1 MATRIX)
C B = NAME DES 2. EINGABEVEKTORS (M-MAL-1 MATRIX)
C R = SKALAR-PRODUKT
C M = ZAHL DER ZEILEN VON A = ZAHL DER ZEILEN VON B
C
C
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(M,1),B(M,1)
C
R = 0.D0
C
DO 10 I = 1,M
10 R = R + A(I,1)*B(I,1)
RETURN
END
c
C
SUBROUTINE ANGLE(V,A,D)
C
C BERECHNUNG VON ALPHA, DELTA AUS DEM VEKTOR DER RICHTUNGSKOSIN.
C
C EINGABE : V = VEKTOR DER RICHTUNGSKOSINUSSE
C AUSGABE : A,D = ALPHA, DELTA IM BOGENMASS
C
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 V(3,1)
C
C
PI = 3.141592653589793D0
C
C
C UEBERGANG ZUM EINHEITSVEKTOR.
C DIES IST DANN NOTWENDIG, WENN VERSEHENTLICH NICHT DIE
C RICHTUNGSKOSINUSSE SELBST, SONDERN EIN VEKTOR, DER EIN
C VIELFACHES DAVON IST, EINGEGEBEN WURDE
C
CALL SCLPRD (V,V,X,3)
X = DSQRT(X)
DO 10 I = 1,3
10 V(I,1) = V(I,1)/X
C
D = DASIN(V(3,1))
A = DATAN2(V(2,1),V(1,1))
IF ( A.LT.0.D0 ) A = A + 2.D0*PI
RETURN
END
c
C
SUBROUTINE PRCPOS (AOLD,DOLD,ANEW,DNEW,PX)
C
C PRAEZESSIONSUEBERTRAGUNG VON OERTERN
C
C EINGABE: AOLD = REKTASZENSION FUER AUSGANGSAEQUINOX
C DOLD = DEKLINATION FUER AUSGANGSAEQUINOX
C PX = PRAEZESSIONSMATRIX
C AUSGABE: ANEW = REKTASZENSION FUER ENDAEQUINOX
C DNEW = DEKLINATION FUER ENDAEQUINOX
C ALLE WINKEL IM BOGENMASS
C
C
IMPLICIT REAL*8 (A-H,O-Z)
REAL*8 PX(3,3),RCOLD(3,1),RCNEW(3,1)
C
C
C BERECHNUNG DER RICHTUNGSKOSINUSSE
C
CALL DIRCOS (AOLD,DOLD,RCOLD)
C
C PRAEZESSIONSUEBERTRAGUNG
C
C
CALL GMPRD (PX,RCOLD,RCNEW,3,3,1)
C
C BERECHNUNG VON ALPHA, DELTA AUS DEN RICHTUNGSKOSINUSSEN
C
CALL ANGLE (RCNEW,ANEW,DNEW)
C
RETURN
END
c
c
SUBROUTINE GMPRD (A,B,R,N,M,L)
C
C MULTIPLIKATION ZWEIER ALLGEMEINER MATRIZEN, R = A.B
C A = NAME DER 1. EINGABEMATRIX
C B = NAME DER 2. EINGABEMATRIX
C R = NAME DES PRODUKTES R = A.B
C N = ZAHL DER ZEILEN VON A = ZAHL DER ZEILEN VON R
C M = ZAHL DER SPALTEN VON A = ZAHL DER ZEILEN VON B
C L = ZAHL DER SPALTEN VON B = ZAHL DER SPALTEN VON R
C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A(N,M),B(M,L),R(N,L)
C
C
DO 20 I = 1,N
DO 20 J = 1,L
R(I,J) = 0.D0
DO 20 K = 1,M
20 R(I,J) = R(I,J) + A(I,K)*B(K,J)
RETURN
END
c
C
SUBROUTINE SMPY (A,C,R,N,M)
C
C MULTIPLIKATION EINER MATRIX MIT EINEM SKALAR
C
C A = NAME DER EINGABEMATRIX
C C = SKALAR
C R = NAME DER AUSGABEMATRIX, R = C.A
C N = ZAHL DER ZEILEN VON A
C M = ZAHL DER SPALTEN VON A
C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A(N,M),R(N,M)
C
C
DO 10 I = 1,N
DO 10 J = 1,M
10 R(I,J) = C*A(I,J)
RETURN
END
c
C
SUBROUTINE GMADD (A,B,R,N,M)
C
C ADDITION ZWEIER MATRIZEN : R = A + B
C
C A = NAME DES 1. SUMMANDEN
C B = NAME DES 2. SUMMANDEN
C R = NAME DER SUMME: R = A + B
C N = ZAHL DER ZEILEN VON A, B, R
C M = ZAHL DER SPALTEN VON A, B, R
C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION A(N,M),B(N,M),R(N,M)
C
C
DO 20 I = 1,N
DO 20 J = 1,M
20 R(I,J) = A(I,J) + B(I,J)
RETURN
END
c
subroutine checklimits (num,delta,xmag,xcol,ier1,ier2,ier3)
c
c Check whether the star's declination and magnitude are
c within the limits of application
c
c Return-Codes
c
c ierr: Error code of the subroutine
c ier1=ier2=ier3=0 -> everything is o.k.
c ier1=1 -> Magnitude is outside the limits for which
c the systematic differecnes can be applied
c ier2=1 -> declination is outside the zone where the
c catalogue comparison can be applied
c ier3=1 -> the star's colour is less than -0.5 or
c greater than 3.0.
c
implicit real*8 (a-h,o-z)
c
common/xlimits/ xdmin,xdmax,xmmin,xmmax
c
c if (1.gt.0) goto 900
c
c
trd = .174532925199433D-01
c
delta2 = delta/trd
c
c Check that the input data are within the availability of
c the catalogue comparison.
c
ier1 = 0
ier2 = 0
ier3 = 0
c
if ( (xmmin .gt.xmag ) .or. (xmag.gt.xmmax ) ) ier1 = 1
if ( (xdmin .gt.delta2) .or. (delta2.gt.xdmax) ) ier2 = 1
if ( (-0.5d0.gt.xcol ) .or. (xcol .gt.3.d0 ) ) ier3 = 1
c
if (ier1.eq.1) write(*,'(
f ''Magnitude not within limits of comparison: '',/,
f ''Star Number'', i8,/,''Magnitude'',f10.2)') num,xmag
c
if (ier2.eq.1) write(*,'(
f ''Declination not within limits of comparison: '',/,
f ''Star Number'', i8,/,''Declination'',f10.2)') num,delta2
c
if ( (ier1.eq.1) .and. (xmag.gt.xmmax) ) write(*,'(
f ''The computation of the systematic differences is made '',/,
f ''with the limiting magnitude '',f8.2,//)') xmmax
if ( (ier1.eq.1) .and. (xmag.lt.xmmin) ) write(*,'(
f ''The computation of the systematic differences is made '',/,
f ''with the limiting magnitude '',f8.2,//)') xmmin
c
if ( ier3.eq.1 ) write(*,'(
f ''The stars colour '',f10.2,'' is outside the allowed '',
f ''interval from -0.5 through 3.0 '',/,
f ''No colour effects are applied for this star '',/)') xcol
if ( ier2.eq.2 ) write(*,'(
f ''The star is skipped'',//)')
c
900 return
end
c
subroutine output(lun,num,xa,xd,xmya,xmyd,epa,epd,xmag)
c
implicit real*8 (a-h,o-z)
c
c The following output data are available (paramter list):
c
c num = Star number in the catalogue
c xa = right ascension in radian
c xd = declination in radian
c xmya = proper motion in right ascension in radian per century
c xmyd = proper motion in declination in radian per century
c xmag = apparent magnitude
c
c sig_sysdif_pos(1) = mean error of the systematic difference in
c right ascension in mas
c
c sig_sysdif_pos(2) = mean error of the systematic difference in
c declination in mas
c
c sig_sysdif_eb(1) = mean error of the systematic difference in
c proper motion in right ascension in mas/yr
c
c sig_sysdif_eb(2) = mean error of the systematic difference in
c proper motion in declination in mas/yr
c
real*8 sig_sysdif_pos(4),sig_sysdif_eb(4)
c
character*12 a12,d12
c
character*1 vorza,vorzd,plus,minus
c
common/sigsys/sig_sysdif_pos,sig_sysdif_eb
c
tra = .261799387799149D+00
trd = .174532925199433D-01
trma = .727220521664304D-04
trmd = .484813681109536D-05
c
plus = '+'
minus = '-'
c
call radcha(xa,'HH MM SS.SSS',a12,ier,6)
call radcha(xd,'VDD MM SS.SS',d12,ier,6)
xmya = xmya/trma
xmyd = xmyd/trmd
c
vorza = plus
if (xmya.lt.0.d0) vorza = minus
xmya = dabs(xmya)
vorzd = plus
if (xmyd.lt.0.d0) vorzd = minus
xmyd = dabs(xmyd)
c write(lun,2300) num,a12,vorza,xmya,d12,vorzd,xmyd,xmag,
c f sig_sysdif_pos(1),sig_sysdif_pos(2),
c f sig_sysdif_eb (1),sig_sysdif_eb (2)
c write(lun,2300) num,a12,vorza,xmya,d12,vorzd,xmyd,xmag
c
c2300 format(i6,2x,a,2x,a,f6.3,2x,a,2x,a,f6.2,f10.2,4f10.2)
c
c
c
c xt = 70.d0
epa = epa - 1900.d0
epd = epd - 1900.d0
c
write(lun,2300) num,a12,vorza,xmya,d12,vorzd,xmyd,epa,epd,xmag
c
c2300 format(i5,1x,a,2x,a,f6.3,2x,a,2x,a,f6.2,f10.2)
2300 format(i6, a,2x,a,f6.3,2x,a,2x,a,f6.2,10x,
f 36x,f7.2,11x,f7.2,13x,f6.2)
c2300 format(i5,1x,a,2x,a,f6.3,2x,a,2x,a,f6.2,f10.2,
c f 36x,f7.2,11x,f7.2)
c
return
end
c
c
subroutine liescat(lun,num,alpha,delta,xmya,xmyd,
F epa,epd,vel ,par,xmag,xcol,ier)
c
c This subroutine reads the catalogue data
c and returns them to the main program.
c Here: the FK5 is read as an example
c
c The user has to adjust this subroutine to his catalogue format.
c
c The following data have to be provided are returned:
c
c alpha : catalogue right ascension (in radians)
c delta : catalogue declination (in radians)
c xmya : proper motion in right ascension (in radians per century)
c xmyd : proper motion in declination (in radians per century)
c epa : mean epoch, for which the right ascension is given (years)
c epd : mean epoch, for which the declination is given (years)
c par : parallaxe of the star (in radians)
c vel : radial velocity
c xmag : The star's apparent magnitude (magnitude in the HIPPARCOS
c system should be given preferentially)
c xcol : The star's colour index B-V in the Johnson system
c
c Units: radians and km/s for the radial velocity
c If no radial velocity or paraalaxe is available
c it should be defined as zero.
c
c
implicit real*8 (a-h,o-z)
c
character*220 feld
character*12 a12,d12
c
c data minus /'-'/
c
c
tra = .261799387799149D+00
trd = .174532925199433D-01
trma = .727220521664304D-04
trmd = .484813681109536D-05
c
ier = +1
c
c Read a star
c
read(lun,'(a)',end=900) feld
c
read(feld,'(i4, 2x,a,f9.3,2x,a,f9.2,46x,
f 2(f7.2,11x),f8.2,10x,f6.5,f9.2)',end=900)
f num,a12,xmya,d12,xmyd,
f epa,epd,xmag1,par,vel
read(feld,'(205x,f8.2)') xmag2
c
xmag = xmag1
xmag = xmag2
xcol = 99. ! No colour is given, no colour effects exis
c
c
call charad(a12,'HH MM SS.SSS',alpha,ier,6)
call charad(d12,'VDD MM SS.SS',delta,ier,6)
c
xmya = xmya*trma
xmyd = xmyd*trmd
par = (par/3600.d0)*trd
c
c Epochs are those, for which the data are really given
c
epa = 2000.d0
epd = 2000.d0
c
return
900 ier = -1
return
end
c
c
Back to the Index