n wenjian tiqu.f90

计算卫星位置的fortran程序

program nwenjian
implicit none
integer prn,nian,yue,ri,i,nj,k,toc,shi,fen,miao,WN,TOW
integer::status=0
REAL*8 wpian,wpiao,wpsu,iode,crs,on,m0,cuc,e,cus,sqrta,toe,cic,omEga,cis,i0,crc,w,omegad
REAL*8 idot,l2,l2p,wjing,wjian,tgd,iodc,fassk,nhqj,x,y,z,gm,GPSZH,T
REAL*8 n,n0,M,F,U,SU,SR,SI,R,IL,JD,EK,WE,L,XK,YK,ZK
character(len=79):: string,forr
character(13):: chara,for
logical::tuichu=.false.
parameter(gm=3.9860047e14,We=.72921151467E-04)    !gm地球引力常数,WE地球自传角速度
for="END OF HEADER"

!    open(unit=10,file="33333.txt",iostat=status)
    open(unit=11,file="f17390700.nav",status="replace",iostat=status)
        do while(.true.)

	       if(.not.tuichu)then     !遇到给定值时开始录如数据
              read(10,'(a79)',iostat=status)string
	              chara=string(61:73)
	   	         if(chara==for)then  !判断条件，符合即开始录入数据
                   tuichu=.true.
		         endif
	       end if
		   if(tuichu)then    

   read(10,'(i2,1x,i2.2,4(1x,i2),f5.1,3d19.12/,3x,4d19.12/,3x,4d19.12/,3x,4d19.12/,&
   3x,4d19.12/,3x,4d19.12/,3x,4d19.12/,3x,2d19.12)',iostat=status)prn,nian,yue,&
   ri,shi,fen,miao,wpian,wpiao,wpsu,iode,crs,on,m0,cuc,e,cus,sqrta,toe,cic,omEga,cis,&
   i0,crc,w,omegad,idot,l2,GPSZH,l2p,wjing,wjian,tgd,iodc,fassk,nhqj
	   IF(STATUS/=0)EXIT
   write(*,'(i2,1x,i2.2,4(1x,i2),f5.1,3d19.12/,3x,4d19.12/,3x,4d19.12/,3x,4d19.12/,&
   3x,4d19.12/,3x,4d19.12/,3x,4d19.12/,3x,2d19.12)')prn,nian,yue,ri,shi,fen,&
   miao,wpian,wpiao,wpsu,iode,crs,on,m0,cuc,e,cus,sqrta,toe,cic,omega,cis,i0,&
   crc,w,omegad,idot,l2,gpszh,l2p,wjing,wjian,tgd,iodc,fassk,nhqj

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    prn=06                              !       试验数据                        ! 
!	 nian=1997                        !  GPS测量原理及应用P33 页                  ! 
!	 yue=11                                                                      ! 
!	 ri=9                                                                        ! 
!	 shi=2                                                                       ! 
!	 fen=0                                                                       
!	 miao=0                                                                            
!	 wpian=-0.231899321079D-06                                                   
!	 wpiao=0                                                                      
!	 wpsu=0
!	 iode=0.970000000000D+02
!	 crs=0.406250000000D+01
!	 on=0.451411660250D-08
!	 M0=-0.290282040486D+00
!	 cuc=0.189989805222D-06
!	 e=0.678421219345D-02
!	 cus=0.912137329578D-05
!	 sqrta=0.515365263176D+04
!	 toe=0.720000000000D+04
!	 cic=0.130385160446D-07
!	 omEga=-0.137835982556D+01
!	 cis=0.949949026108D-07
!	 i0=0.958512160302D+00
!	 crc=0.201875000000D+03
!	 w=-0.258419417299D+01
!	 omegad=-0.819426989566D-08
!	 idot=-0.253939149013D-09
!	 Gpszh=0.931000000000D+03
!	 wjing=0.700000000000D+03
!	 wjian=0
!	 tgd=0.186264514923D-08
!	 iodc=0.353000000000D+03

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


!       WRITE(*,'(A,D18.12)')"M0=",M0
!第一轨道参数   平均角速度        
	   N0=SQRT(GM)/(SQRTA**3)
	   N=N0+ON
	   WRITE(*,'(A,I2,A,D18.12)')"第",PRN,"颗卫星的平均角速度n：",n

!计算卫星的平近点角
!       IF(NIAN>=80)THEN   !大于80年，肯定是1980年之后，年+1900，
!	     NIAN=NIAN+1900
!	   ELSE 
!	     NIAN=NIAN+2000        !小于80，则年肯定是2000年之后，所以年+2000
!	   END IF
       CALL G_GPS(NIAN,YUE,RI,SHI,FEN,MIAO,WN,T)
!	     IF(WN/=GPSZH)THEN
!	        WN=WN-1
!		    T=T+604800
!	     END IF
!	   WRITE(10,'(D18.12)')M0
       DO I=0,39
	      T=T+180*I
       M=M0+N*(T-TOE)
	          WRITE(*,'(A,D18.12)')"平近点角M：",M
!计算偏近点角
       CALL P_JDJ(M,E,EK)
	   WRITE(*,'(A,D18.12)')"偏近点角E：",EK

!计算真近点角
       F=ATAN(SQRT(1-E*E)*SIN(EK)/(COS(EK)-E))
       WRITE(*,'(A,D18.12)')"真近点角F：",F

!计算生交距角U'
       U=W+F
       WRITE(*,'(A,D18.12)')"生交距角U：",u

!计算摄动改正项SU,SR,SI
       SU=CUC*COS(2*U)+CUS*SIN(2*U)
	   SR=CRC*COS(2*U)+CRS*SIN(2*U)
	   SI=CIC*COS(2*U)+CIS*SIN(2*U)
	   WRITE(*,'(3(A,D21.15))')"SU:",SU,"SR:",SR,"SI：",SI
	   	   !对u',r',i0'进行摄动改正
	       U=U+SU
           R=SQRTA*SQRTA*(1-E*COS(EK))+SR
		   IL=I0+SI+OMEGAD*(T-TOE)      !IL即i0
       WRITE(*,'(3(A,D21.15))')"摄动改正项U':",U,"R:",R,"I：",IL

!计算卫星在轨道面坐标系中的位置:
       XK=R*COS(U)
	   YK=R*SIN(U)
       WRITE(*,'(A/,A,D20.12,2X,A,D20.11)')"卫星在轨道平面直角坐标系中的坐标为：","X=",XK,"Y=",YK

!计算观测瞬间升交点的经度L:
       L=OMEGA+(OMEGAD-WE)*T-OMEGAD*TOE
       WRITE(*,'(A,D19.12)')"升交点经度L:",L

!计算卫星在瞬时地球坐标系中的位置
       X=XK*COS(L)-YK*COS(IL)*SIN(L)
	   Y=XK*SIN(L)+YK*COS(IL)*COS(L)
	   Z=YK*SIN(IL)
       WRITE(11,'(3(D18.12,2X))')X,Y,Z
	   END DO

       WRITE(*,*)"****************************************"
!	   IF(STATUS/=0)EXIT	  
!      ENDIF 
!    end do
!  close(10)
  close(11)
 stop
end