! ********************************************************* ! mainpro vers. 9 ! May 2000 ! this program computes at once mean and proper elements ! version with OEF (OrbFit compatible) input ! ********************************************************* PROGRAM prop9 USE massmod USE short9 ! masses IMPLICIT NONE ! max no. of entries INTEGER, PARAMETER :: numx=20000 ! ============================================================= ! el,em: 1=semimaj. ax. 2= ecc. 3=inclin. ! 4=node 5=arg. peri 6=mean anomaly ! planetary elements; second index is 1=jup 2=sat DOUBLE PRECISION pel(6,2) ! asteroid keplerian elements (osculating, mean) DOUBLE PRECISION el(6),em(6) ! end of file logical eof ! file names ! these are for input planets and baric character*72 filpla,filbar ! real input/output file names character*72 filast,filfil,filpro INTEGER len ! options INTEGER inp,nmax ! input to short8s subroutine DOUBLE PRECISION tt CHARACTER*9 nam0 ! reference time for long9 DOUBLE PRECISION t0 ! outputs from short8 subroutine INTEGER iqcm,it,ierr DOUBLE PRECISION delta ! osculating elements quality code, magnitude just copied INTEGER iqco DOUBLE PRECISION hm ! counters INTEGER kr,krdon ! output proper elements:frequencies, actions, angles DOUBLE PRECISION gp,sp,g0,s0,pre,prsini,apip,atep,apiv,atev INTEGER ngper,ngnod ! quality codes INTEGER iqce,iqcf,irfl,irfl2,ipas ! unit conversion DOUBLE PRECISION conv ! unit numbers INTEGER iun1,iun2,iun3,iun9,iun10,iun11,iun19,iun21,iun18,iun8 & & ,iun7,iun17, iundone ! ******************************************************** ! use orbfit library CALL libini ! choice of input format, options: CALL filopn(iun1,'prop.opt','old') call skip (iun1,1) ! inp=1 from catalogue; inp=2 numerical integration; inp=3 filtered dat call reaint(iun1,'inp',inp) if(inp.le.0.or.inp.gt.3)then write(*,*)' unknown input option, inp= ',inp stop endif ! control on number of asteroids call reaint(iun1,'nmax',nmax) ! ---------------------------------------------------------------------- ! input/output file names ! input file: asteroid osculating elements call reastr(iun1,'filast',filast) ! input file for mean (filtered) elements: ! (only for inp=3, otherwise dummy) call reastr(iun1,'filfil',filfil) ! output file: beginning of name call reastr(iun1,'filpro',filpro) ! ---------------------------------------------------------------------- ! planetary data for short periodic perturbations and barycentric correc call skip(iun1,1) ! input file: planets call reastr(iun1,'filpla',filpla) ! input file: barycenter of the inner solar system ! ibar=1 barycentric correction, ibar=0 no correction call reastr(iun1,'filbar',filbar) ! end input options close(iun1) ! ********************************************************************* ! initialisations, file opening, for short and long ! ********************************************************************** call openfi(inp,t0,filast,filfil,filpro,filpla,filbar, & & iun2,iun3,iun9,iun10,iun11,iun19,iun21,iun18,iun7,iun17) ! write(*,*)iun8,iun18,iun7,iun17 ! ********************************************************* ! loop on input asteroid data begins here ! counter for mean/proper elements krdon=0 DO 10 kr=1,nmax ! input osculating elements if(inp.eq.1.or.inp.eq.2)then call iosho(kr,inp,iun8,iun7,iun17,nam0,t0,tt,el, & & pel,eof,iqco,hm,iun19) if(kr.eq.100*(kr/100))then if(inp.eq.1)then write(*,*)nam0,' read no ',kr elseif(inp.eq.2)then write(*,*)tt,' read no ',kr endif endif if(eof)goto 99 ! ********************************************************* ! excluding the non-main-belt asteroids. if(inp.eq.1)then if ((el(1)*(1.d0-el(2))).lt.1.1d0.or.el(1).gt.3.8d0 & & .or.el(1).lt.1.8d0)then goto 10 endif endif ! counter for mean elements krdon=krdon+1 ! *************************************************** CALL short(el,pel,nam0,tt,inp,em,iqcm,it,delta,ierr,iun19,iun21) ! ******************************************************** ! output of the mean elements if(inp.eq.1)then write(iun21,911)nam0,it,delta 911 format(a9,i3,d12.4) write (iun19,410) nam0,em(6),em(5),em(4),em(3),em(2), & & em(1),iqcm,ierr,iqco,hm 410 format(a9,4f10.5,f11.7,f11.7,i3,i3,i3,1x,f5.2) elseif(inp.eq.2)then write(iun21,912)tt,it,delta 912 format(f12.2,i3,d12.4) write (iun19,411)tt,em(6),em(5),em(4),em(3),em(2),em(1), & & iqcm,ierr 411 format(f12.2,4f10.5,f10.7,f12.7,i3,i3) endif ! ******************************************************** else ! filtered mean elements to be read for inp=3 ! input from output of orbit8v, filtered ! only one orbit in the input, equinoctal call inpfil(iun18,tt,em,eof) ! iqcm is not given, it is assumed that filtering has ! eliminated all the short periods anyway iqcm=0 ! end job and counters if(eof)goto 99 if(kr.eq.(kr/100)*100)then write(*,*)tt,' read ',kr endif ! ******************************************************** ! mean elements are available endif ! compute proper elements call long9(em,nam0,tt,t0,inp, & & gp,sp,pre,prsini,iqce, & & apip,atep,g0,s0,ipas,iqcf, & & irfl,irfl2, & & ngper,apiv,ngnod,atev,conv, & & iun2,iun3,iun9,iun10,iun11) ! output proper elements call selout(inp,tt,nam0, & & gp,sp,em(1),pre,prsini,iqce,iqcm,iqco, & & hm,apip,atep,g0,s0,ipas,iqcf, & & irfl,irfl2, & & ngper,apiv,ngnod,atev,conv,iun2,iun3,iun9,iun11) ! counter if(krdon.eq.100*(krdon/100))then if(inp.eq.1)then write(*,*)nam0,' mean el. done no ',krdon elseif(inp.eq.2.or.inp.eq.3)then write(*,*)tt,' mean el. done no ',krdon endif endif ! *************************************************** ! end loop on asteroid catalog 10 ENDDO ! displaying the information on job termination. 99 continue write(*,*)krdon,' mean and proper elements generated ' ! proper elements output files (always) call filclo(iun2,' ') call filclo(iun9,' ') call filclo(iun10,' ') ! discarded proper elements (only for catalog input) IF(inp.eq.1)call filclo(iun3,' ') IF(inp.eq.1)call filclo(iun11,' ') ! mean elements output files (if mean elements are computed) IF(inp.eq.1.or.inp.eq.2)call filclo(iun19,' ') IF(inp.eq.1.or.inp.eq.2)call filclo(iun21,' ') ! input planetary data (if mean elements are computed) IF(inp.eq.1.or.inp.eq.2)call filclo(iun17,' ') IF(inp.eq.1.or.inp.eq.2)call filclo(iun7,' ') ! input filtered data (if mean elements are not computed) IF(inp.eq.3)call filclo(iun18,' ') CALL filopn(iundone,'prop9.done','unknown') CALL filclo(iundone,' ') END PROGRAM prop9