trans.for

坐标转换的实现 可批量进行高精度转换,实现了测量数据的换带计算及高程投影的进行.

        PROGRAM MAIN
        real*8 CMQD,CMHD,CMQ,CMH,HQ,HH,ZXQ,ZYQ,ZXH,ZYH
        real*8 X,Y,B,L,XXX,YYY,XX,YY,bbb,lll,BB,LL
	integer NUM
        character*13 INPUT,OUTPUT,NAME
        character*50 XMNAME
        WRITE(*,*) (' ')
        WRITE(*,*) ('     ==================================')
        WRITE(*,*) ('          PROJECTION TRANSFORMATION')
        WRITE(*,*) ('                    V1.0')
        WRITE(*,*) ('                HU KAICHENG')
        WRITE(*,*) ('     ==================================')
        WRITE(*,*) (' ')
        WRITE(*,*) ('INPUT FILENAME')
*        READ(*,10)INPUT
        WRITE(*,*) ('OUTPUT FILENAME')
*        READ(*,10)OUTPUT
        OPEN(1,FILE='IN.TXT')
        OPEN(2,FILE='OUT.TXT')
        OPEN(3,FILE='TRANS.TXT')
        READ(3,*)XMNAME
        READ(3,*)CMQD,CMHD
        READ(3,*)HQ,HH
        READ(3,*)ZXQ,ZYQ,ZXH,ZYH
        call DFMTOD(CMQD,CMQ)
        call DFMTOD(CMHD,CMH)
        WRITE(2,30) ('项  目  名  称:'),XMNAME
 30     FORMAT(10X,A15,A50)
        WRITE(2,31) ('中 央 子 午 线:'),cmhd
 31     FORMAT(10X,A15,F10.6)
        WRITE(2,32) ('投影抵偿高程面:'),HH,('m') 
 32     FORMAT(10X,A15,F9.3,A1)  
        call XYTOBL(ZXH,ZYH-5D5,BB,LL) 
        LL=LL+CMH
        CALL DTODFM(BB,BBB)
        CALL DTODFM(LL,LLL)
        WRITE(2,33) ('投  影  中  心:'),bbb,lll
 33     FORMAT(10X,A15,f8.4,2X,f8.4)
        WRITE(2,*) (' ')
        READ(1,*) NUM
        do 1000 i=1,num
        READ(1,*) NAME,XXX,YYY
        WRITE(*,*)NAME,XXX,YYY
        call TOUYING(HQ,ZXQ,ZYQ,XXX,YYY,0D0,X,Y)
        call XYTOBL(X,Y-5D5,B,L)
        L=L+CMQ
        call BLTOXY(CMH,B,L,X,Y)
        Y=Y+5D5
        call TOUYING(0D0,ZXH,ZYH,X,Y,HH,XX,YY)
        WRITE(2,20)name,XX,YY 
 1000   continue
 10	FORMAT(A12)
 20	FORMAT(1X,A13,1X,F11.3,3X,F10.3)
        END


        SUBROUTINE XYTOBL(X,Y,B,L)
        real*8 X,Y,B,L,PI,E2,C
        real*8 BF,F1,F2,F3,F4,F5,F6,F7
        real*8 TF,NF2,N
        real*8 B1,B2,B3,l1,l2
        E2=0.0067385254147D0
        C=6399698.90178271D0
        PI=3.14159265358979323846D0
        F1=27.11115372595D0
        F2=9.02468257083D0
        F3=0.00579740442D0
        F4=0.00043532572D0
        F5=0.00004857285D0
        F6=0.00000215727D0
        F7=0.00000019404D0
        BF=F1+F2*(X/1.0D6-3)-F3*(X/1.0D6-3)**2
        BF=BF-F4*(X/1.0D6-3)**3+F5*(X/1.0D6-3)**4
        BF=BF+F6*(X/1.0D6-3)**5-F7*(X/1.0D6-3)**6
        TF=DTAN(BF*PI/1.8D2)
        NF2=E2*(DCOS(BF*PI/1.8D2))**2
        N=Y*SQRT(1+NF2)/C
        B1=5.0D0+3.0D0*TF**2+NF2-9.0D0*NF2*TF**2
        B2=0.25D0*(6.1D1+9.0D0*TF**2+4.5D0*TF**4)*N**6
        B3=9.0D1*N*N-7.5D0*B1*N**4+B2
        B=BF-(1+NF2)*TF*B3/PI
        L1=3.0D1*(1.0D0+2.0D0*TF**2+NF2)*N**3
        L2=1.5D0*(5.0D0+2.8D1*TF**2+2.4D1*TF**4)*N**5
        L=(1.8D2*N-L1+L2)/(PI*COS(BF*PI/1.8D2))
        END

        SUBROUTINE BLTOXY(CL,B,L,X,Y)	
        real*8 B,L,CL,PI,E2
	real*8 XXX,LL,SIB,COB
	real*8 T,N2,N,C,M
	real*8 C0,C1,C2,C3,C4
	real*8 X,Y,XX1,YY1,XX2,YY2,XX0
        PI=3.14159265358979323846D0
        E2=0.0067385254147D0
        C=6399698.90178271D0
        c0=111134.8611D0
        c1=32005.7799D0
        c2=133.9238D0
        c3=0.6976D0
        C4=0.0039D0
        LL=L-CL
        SIB=DSIN(B*PI/180.0D0)
        COB=DCOS(B*PI/180.0D0)
        T=DTAN(B*PI/180.0D0)
        N2=E2*COB**2
        N=C/(SQRT(1+N2))
        M=COB*LL*PI/180.0D0
        XXX=C0*B-(C1*SIB+C2*SIB**3+C3*SIB**5+C4*SIB**7)*COB
        XX0=5.0D0-T*T+9.0D0*N2+4.0D0*N2*N2
        XX1=M**2/2.0+(XX0*M**4.0)/24.0D0
        XX2=((61.0D0-58.0D0*T*T+T**4)*M**6)/720.0D0
        X=XXX+N*T*(XX1+XX2)
        YY1=1.0D0-T*T+N2
        YY2=5.0D0-18.0D0*T*T+T**4+14.0D0*N2-58.0D0*N2*T*T
        Y=N*(M+YY1*M**3/6.0D0+YY2*M**5/120.0D0)
        END
        

        SUBROUTINE DFMTOD(A,B)
        REAL*8 A,b,K1,K2,K3
        A=A+1D-14
        K1=DINT(A)
        K2=DINT(A*100D0)-DINT(A)*1D2
        K3=A*10000D0-DINT(A*100D0)*1D2
        B=K1+K2/6D1+K3/3.6D3
        END

        SUBROUTINE DTODFM(U,V)
        REAL*8 U,V,KK1,KK2,KK3
        KK1=DINT(U)
        KK2=DINT((U-KK1)*6D1)
        KK3=(U-KK1-KK2/6D1)*3.6D3
        V=KK1+KK2/1D2+KK3/1D4
        END
        
        SUBROUTINE TOUYING(HQ,XZX,YZX,XQ,YQ,H,XH,YH) 
        REAL*8 HQ,XZX,YZX,XQ,YQ,H,XH,YH,RE
        RE=6.371D6
          XH=Xzx+(XQ-XZX)*((RE+H)/(RE+HQ))
          YH=Yzx+(YQ-YZX)*((RE+H)/(RE+HQ))   
        END