! ================================================================== ! input/output routines for longit8x; version 8 ! ================================================================== ! outpro ! print proper elements, frequencies ! ================================================================== ! in this version, proper elements are printed anyway; ! however, the quality code contains warnings and some are in a ! separate file; also added quality code for mean elements SUBROUTINE outpro(inp,t,nam0,aa,pre,prsini,gsy,ssy,iqce, & & api1,ate1,g0sy,s0sy,iqcf,npas,ngper,apiv,ngnod,atev & & ,ipro,iang,iqcm,iqco,hm,irfl,irfl2) USE fund_const IMPLICIT NONE INTEGER, INTENT(IN) :: inp ! input mode ! string for alphanumeric asteroid name CHARACTER*9, INTENT(IN) :: nam0 DOUBLE PRECISION, INTENT(INOUT ) :: aa,pre,prsini ! proper elems. DOUBLE PRECISION, INTENT(INOUT) :: api1,ate1,apiv,atev,g0sy,s0sy,gsy,ssy ! proper angles, frequencies DOUBLE PRECISION, INTENT(IN) :: hm,t ! magnitude, epoch INTEGER, INTENT(IN) :: iqcm,iqco,irfl,irfl2 ! quality codes, flags INTEGER, INTENT(INOUT) :: iqce,iqcf INTEGER, INTENT(IN) :: npas,ngper,ngnod ! no. iter, no. revolutions INTEGER, INTENT(IN) :: ipro,iang ! output units ! END INTERFACES INTEGER ir DOUBLE PRECISION api, ate ! safeties against format errors (absurd cases) if(pre.gt.0.9999d0)pre=0.9999d0 if(prsini.gt.0.9999d0)prsini=0.9999d0 if(gsy.gt.9999.d0)gsy=9999.d0 if(gsy.lt.-9999.d0)gsy=-9999.d0 if(ssy.gt.9999.d0)ssy=9999.d0 if(ssy.lt.-9999.d0)ssy=-9999.d0 if(g0sy.gt.9999.d0)g0sy=9999.d0 if(g0sy.lt.-9999.d0)g0sy=-9999.d0 if(s0sy.gt.9999.d0)s0sy=9999.d0 if(s0sy.lt.-9999.d0)s0sy=-9999.d0 if(inp.eq.1.or.inp.eq.4)then ! **************************************************** ! catalogue of proper elements if(iqce.lt.0)iqce=0 if(iqcf.lt.0)iqcf=0 if(irfl.eq.0)then if(iqce.gt.10)then ir=90+(iqce-10)/3 else ir=0 endif else ir=irfl endif api=api1*360/dpig if(api.gt.360)api=api-360.d0 ate=ate1*360/dpig if(ate.gt.360)ate=ate-360.d0 write (ipro,200) nam0,aa,pre,prsini,gsy,ssy, & & ir,iqcm,iqco,hm 200 format (a9,f8.5,1x,f5.4,1x,f5.4,2f10.3,i4,i3,i3,1x,f5.2) write(iang,400)nam0,api,ate & & ,g0sy,s0sy,npas,iqcf,iqce,irfl2 400 format(a9,2f9.3,2f10.3,4i4) elseif(inp.eq.2.or.inp.eq.3)then ! **************************************************** ! output of numerical integration ! the proper angles are computed for current time, with number of ! revolutions accounted for; output in degrees call prngup(ate1,atev,ngnod) ate1=ate1*360/dpig ate1=ate1+ngnod*360 ! for reduction to t=t0: -(t-t0)*ssy/3600 call prngup(api1,apiv,ngper) api1=api1*360/dpig api1=api1+ngper*360 if(irfl.eq.0)then if(iqce.gt.10)then ir=90+(iqce-10)/3 else ir=0 endif else ir=irfl endif ! for reduction to t=t0: -(t-t0)*gsy/3600 write (ipro,201) t,aa,pre,prsini,gsy,ssy,ir,iqcm 201 format (f12.2,f12.8,2f9.6,2f12.4,i4,i3) write (iang,202) t,api1,ate1,g0sy,s0sy,npas,iqcf, & & iqce,irfl2 202 format (f12.2,2f15.5,2f12.4,4i4) endif END SUBROUTINE outpro ! --------------------------------------------------------------- ! input of filtered output and conversion to keplerian SUBROUTINE inpfil(iun,tt,em,eof) USE fund_const IMPLICIT NONE INTEGER iun DOUBLE PRECISION em(6),tt LOGICAL eof DOUBLE PRECISION eh,ek,eq,ep,awt read(iun,*,end=99)tt,em(1),eh,ek,ep,eq eof=.false. if(eh.eq.0.d0.and.ek.eq.0.d0)then awt=0.d0 else awt=atan2(eh,ek)*degrad endif if(ep.eq.0.d0.and.eq.eq.0.d0)then em(4)=0.d0 else em(4)=atan2(ep,eq)*degrad endif em(5)=awt-em(4) em(2)=sqrt(eh**2+ek**2) em(3)=atan(sqrt(ep**2+eq**2))*2.d0*degrad RETURN 99 eof=.true. END SUBROUTINE inpfil ! ================================================================== ! openfi ! open and intest routine for longit7 ! ********************************************************************* SUBROUTINE openfi(inp,t0,filast,filfil,filpro,filpla,filbar & & ,iun2,iun3,iun9,iun10,iun11,iun19,iun21,iun18,iun7,iun17) IMPLICIT NONE ! input type, file names INTEGER inp CHARACTER*60 filpro,filast,filfil,filpla,filbar ! output ! reference time, only for inp=3 DOUBLE PRECISION t0 ! unit numbers INTEGER iun2,iun3,iun9,iun10,iun11,iun19,iun21,iun18,iun7,iun17 ! end interface INTEGER len,nast,ilce CHARACTER*60 filang,filres,filmea,filnam,fildis,filadi,filrem,filprp CHARACTER*3 coox,unitx,sysx CHARACTER*6 refx CHARACTER*4 nomast(20) CHARACTER*60 comast CHARACTER*100 colhea INTEGER ibar ! open input files IF(inp.eq.1.or.inp.eq.2)THEN ! input either catalog or time series, but of osculating elements anyway ! open files for input of osculating elements call rmbl(filast,len) IF(inp.eq.1)THEN CALL oporbf(filast(1:len), -1) ! ELSEIF(inp.eq.2)THEN ! call filopn(iun8,filast(1:len),'old') ENDIF ! open input file for planetary data call rmbl(filpla,len) call filopn(iun7,filpla,'old') ! open input file for barycenter of inner solar system if(inp.eq.2)then ibar=0 write(*,*)' no barycentric correction in time series' elseif(inp.eq.1)then ibar=1 write(*,*)' baricenter from file ',filbar call rmbl(filbar,len) call filopn(iun17,filbar,'old') elseif(inp.eq.3)then ibar=22 else write(*,*)' unknown input type, inp=',inp stop endif ELSEIF(inp.eq.3)THEN ! input directly from orbit9 (digitally filtered) call rmbl(filfil,len) call filopn(iun18,filfil(1:len),'old') call reahea(iun18,nomast,comast,colhea,coox,sysx,refx, & & unitx,nast,ilce,t0) if(nast.gt.1)then write(*,*)nast,' orbts at once not allowed; use conv8v' stop endif if(coox.ne.'EQU'.or.sysx.ne.'HEL'.or.refx.ne.'INVL1B')then write(*,*)' wrong input coordinates/ref.system, use conv8v' stop endif ELSE WRITE(*,*)' inp =',inp,' not allowed' STOP ENDIF ! open output files IF(inp.eq.1.or.inp.eq.2)THEN ! open files for output of mean elements call rmbl(filpro,len) filmea=filpro(1:len)//'.mea' call filopn(iun19,filmea(1:len+4),'unknown') filrem=filpro(1:len)//'.rem' call filopn(iun21,filrem(1:len+4),'unknown') ENDIF ! open files with names from prop.opt ! these always exist call rmbl(filpro,len) filprp=filpro(1:len)//'.pro' call filopn(iun2,filprp(1:len+4),'unknown') filang=filpro(1:len)//'.ang' call filopn(iun9,filang(1:len+4),'unknown') filres=filpro(1:len)//'.res' call filopn(iun10,filres(1:len+4),'unknown') IF(inp.eq.1)THEN ! discarded cases in a separate file fildis=filpro(1:len)//'.dis' call filopn(iun3,fildis(1:len+4),'unknown') filadi=filpro(1:len)//'.adi' call filopn(iun11,filadi(1:len+4),'unknown') ENDIF ! ********************************************* ! output headers for all input types write(iun2,592) 592 format(' PROPER ELEMENTS MILANI and KNEZEVIC') write(iun9,592) write(iun10,592) if(inp.eq.1)then ! original input from catalog of osculating elements ! no input to long periodic part ! discarded cases in a separate file write(iun3,592) write(iun11,592) ! else ! write(*,*)' input type ',inp,' not supported' ! stop endif END SUBROUTINE openfi ! ================================================================== ! selout ! conversions, output selection, output SUBROUTINE selout(inp,t,nam0, & & gv,sv,aa,pre,prsini,iqce,iqcm,iqco, & & hm,api,ate,g0,s0,ipas,iqcf,irfl,irfl2, & & ngper,apiv,ngnod,atev,conv,iun2,iun3,iun9,iun11) IMPLICIT NONE DOUBLE PRECISION gv,sv,aa,pre,prsini,hm,api,ate, g0,s0,apiv,atev,conv,t INTEGER iqce,iqcm,iqco,iun2,iun3,iun9,iun11,ipas,iqcf,irfl,irfl2, ngper,inp,ngnod LOGICAL discard CHARACTER*9 nam0 ! END INTERFACE DOUBLE PRECISION gsy, ssy, g0sy, s0sy INTEGER ipro, iang ! conversions to arcsec/yr gsy=gv*conv ssy=sv*conv g0sy=g0*conv s0sy=s0*conv call select(inp,aa,pre,prsini,gsy,ssy,iqce,iqcf,iqcm,discard) IF(discard)THEN ipro=iun3 iang=iun11 ELSE ipro=iun2 iang=iun9 ENDIF call outpro(inp,t,nam0,aa,pre,prsini,gsy,ssy,iqce, & & api,ate,g0sy,s0sy,iqcf,ipas,ngper,apiv,ngnod,atev & & ,ipro,iang,iqcm,iqco,hm,irfl,irfl2) END SUBROUTINE selout ! ================================================================== ! elems8 ! if it is meaningful, compute proper e and sin I ! if blow up occurs (e>1 and/or sin I>1) write strong warning ! input: rho(4) proper Poincare' elements ! output: ami,ani proper actions*2; api,ate proper angles ! pre,prsini proper keplerian el ! idisc discarded elements code (see below) SUBROUTINE elems8(rho,ani,ami,api,ate,pre,prsini,idisc) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: rho(4) DOUBLE PRECISION, INTENT(OUT) :: ani,ami,api,ate,pre,prsini INTEGER, INTENT(OUT) :: idisc ! END INTERFACE DOUBLE PRECISION cani, cami ! proper actions (doubled) ani=rho(1)**2+rho(2)**2 ami=rho(3)**2+rho(4)**2 ! proper angles api=datan2(rho(2),rho(1)) ate=datan2(rho(4),rho(3)) ! proper e if(ani.lt.2.d0)then cani=ani-ani**2/4.d0 pre=sqrt(cani) idisc=0 else ! e>1: idisc=1 idisc=1 endif ! proper sin I if(ami/sqrt(1-cani).lt.2.d0)then cami=ami/sqrt(1-cani)-ami**2/4.d0/(1-cani) prsini=sqrt(cami) else ! sin I >1: idisc=2; sin I >1.and.e>1: idisc=3 idisc=idisc+2 endif END SUBROUTINE elems8