三參數坐標轉換 Helmert Transformation
三參數坐標轉換
多年前寫的坐標轉換程式，上週同學告知有bug，當年是以Lahey Fortran編寫程式碼，沒想到Lahey公司2022年底宣布關門，頓時不知如何面對江東父老。
高階語言最麻煩的是處理Unicode，如果編譯器不支援就只能以英文I/O。想破頭兩星期終於找到出路。
PROGRAM HELMERT3B
! =========================================================
! = Modify 2000.12.30 固定輸出入檔案之檔名 =
! = 2001.5.24 修改計算模式 =
! = 90.05.29. 修改轉換參數之輸出有效位數(蕭欽霖之貢獻) =
! = 90.10.14. 輸出轉換參數值, 以供AutoCAD程式使用 =
! = 91.05.09. 修改不符值輸出順序錯誤 =
! = 2004.12.15 compiled by Lahey Fortran 95 for Windows XP =
! = 2014.01.28 將旋轉中心(0,0)移至測區中央 =
! = 適用於兩坐標系統均為二度分帶系統 =
! = 2023.02.19 修正度分秒副程式之bug =
! = 以GNU Fortran Compile =
! ========================================================
IMPLICIT REAL(8)(A-H,O-Z)
REAL(8) XF,YF,XI,YI
CHARACTER INFILE*60,OUTFILE*60,CORFILE*60,JOB*70,PARFILE*60
CHARACTER*10 NAME,PNAME
INTEGER idt(8)
CHARACTER*10 date,time,zone
COMMON AVEX1,AVEY1,AVEX2,AVEY2
! CALL GETDAT(IYEAR,IMONTH,IDAY)
! CALL GETTIM(IHOUR,IMIN,ISEC,I100H)
CALL date_and_time(date,time,zone,idt)
WRITE(*,1000)
1000 FORMAT(20X,'***************************************'/&
20X,'* HELMERT 三參數相似轉換 *'/&
20X,'* 形心對形心坐標轉換 *'/&
20X,'* compiler GNU Fortran 95 *'/&
20X,'* program by ChinaTom *'/&
20X,'* Version 3.0 2023.5. *'/&
20X,'***************************************')
WRITE(*,1002)
1002 FORMAT(/10X,'本程式適用於兩坐標系統均為二度分帶系統,且將旋轉中心'&
'(0,0)移至測區中央')
WRITE(*,1010)
1010 FORMAT(/10X,'輸入資料檔名 <<主檔名即可>>： '$)
READ(*,1020) INFILE
1020 FORMAT(A60)
NN=LEN_TRIM(INFILE)
ICHK=INDEX(INFILE,'.')
IF(ICHK.NE.0) THEN
WRITE(*,1030)
1030 FORMAT(/10X,'<< 輸入錯誤 : 檔名中含附加檔名 >> 程式中止!')
! PAUSE ' 程式中止! 按 ENTER 鍵離開'
STOP
END IF
OUTFILE=INFILE
CORFILE=INFILE
PARFILE=INFILE
INFILE(NN+1:NN+4)= '.DAT'
OUTFILE(NN+1:NN+4)='.DOC'
CORFILE(NN+1:NN+4)='.COD'
PARFILE(NN+1:NN+4)='.PA3'
OPEN(1,FILE=INFILE,STATUS='OLD',IOSTAT=ICHK)
! IF(ICHK.EQ.6416) THEN
IF(ICHK.GT.0) THEN
WRITE(*,1040) INFILE
1040 FORMAT(/10X,'<<錯誤>> : 找不到資料檔 ',A,' 程式中止!')
! PAUSE ' 程式中止! 按 ENTER 鍵離開'
STOP
END IF
OPEN(2,FILE=OUTFILE,STATUS='UNKNOWN')
OPEN(3,FILE=CORFILE,STATUS='UNKNOWN')
OPEN(4,FILE='AA.TMP',STATUS='UNKNOWN',FORM='UNFORMATTED')
OPEN(5,FILE='BB.TMP',STATUS='UNKNOWN',FORM='UNFORMATTED')
OPEN(8,FILE=PARFILE,STATUS='UNKNOWN')
WRITE(2,1000)
NUM=1
NUMRJ=0
READ(1,'(A70)') JOB
10 READ(1,1050,END=20) IC,NAME,YI,XI,YF,XF
1050 FORMAT(I5,A10,4F15.0)
! =============================================================
! = IC : >0 POINT NOT USED, =-99 END OF CONTROL POINTS =
! = NAME : POINT NAME =
! = XI,YI : (X,Y) OF ORIGINAL COORDS. SYSTEM =
! = XF,YF : (X,Y) OF DESTINATION COORDS. SYSTEM =
! =============================================================
IF(IC.GT.0) THEN
WRITE(5) NAME,XI,YI,XF,YF
NUMRJ=NUMRJ+1
GOTO 10
END IF
IF(IC.EQ.-99) GOTO 20
WRITE(4) NAME,XI,YI,XF,YF
NUM=NUM+1
GOTO 10
20 NUM=NUM-1
IF(NUM.LE.1) THEN
WRITE(*,1060)
1060 FORMAT(//10X,'錯誤：控制點點數少於 2點 !')
STOP
END IF
WRITE(2,1070) TRIM(JOB)
1070 FORMAT(/5X,'資料名稱 : ',A)
WRITE(2,1080) TRIM(INFILE)
1080 FORMAT( 5X,'檔案名稱 : ',A)
! WRITE(2,1090) IYEAR,IMONTH,IDAY,IHOUR,IMIN
!1090 FORMAT( 5X,'計算日期時間 : ',I4,'/',I2,'/',I2,2X,I2,':',I2)
WRITE(2,1090) IDT(1),IDT(2),IDT(3),IDT(5),IDT(6),IDT(7)
1090 FORMAT(5X,'計算日期 : ',I4,'/',I2,'/',I2,4X,'時間: ',I2,':',I2,':',I2)
IF(NUMRJ.NE.0) THEN
REWIND 5
WRITE(2,1100)
1100 FORMAT(/1X,'[未加入轉換之已知坐標點位]'/&
11X,'點 名',11X,'轉換前(N,E)坐標',12X,'轉換後(N,E)坐標'/&
3X,72(1H=))
DO 22 I=1,NUMRJ
READ(5) NAME,XI,YI,XF,YF
WRITE(2,1110) I,NAME,YI,XI,YF,XF
1110 FORMAT(4X,I3,3X,A10,2(3X,F11.3,2X,F11.3))
22 CONTINUE
END IF
CALL PARA3(X04,Y04,ANG4,NUM)
! IF(EOF(1)) GOTO 40
WRITE(*,1120)
WRITE(2,1120)
1120 FORMAT(//1X,'[未 知 點 轉 換 成 果]'/&
11X,'點 名',11X,'轉換前(N,E)坐標',12X,'轉換後(N,E)坐標'/&
8X,67(1H=))
30 READ(1,1130,END=40) IC,PNAME,Y,X
1130 FORMAT(I5,A10,2F15.0)
IF(IC.GT.0) GOTO 30
CALL TR3(X,Y,X04,Y04,ANG4,XX,YY)
WRITE(*,1140) PNAME,Y,X,YY,XX
WRITE(2,1140) PNAME,Y,X,YY,XX
1140 FORMAT(10X,A10,2(3X,F11.3,2X,F11.3))
WRITE(3,1150) PNAME,YY,XX
1150 FORMAT(A10,2F15.3)
GOTO 30
40 CONTINUE
CLOSE(4,STATUS='DELETE')
CLOSE(5,STATUS='DELETE')
CLOSE(8)
END
SUBROUTINE PARA3(X0,Y0,ANG,NUM)
! ===========================================================
! = PROGRAM : 3-PARAMETER TRANSFORMATION =
! = (X,Y)-SHIFT, ROTATION =
! =---------------------------------------------------------=
! = X2=X1*COS(ang)-Y1*SIN(ang)+X0 =
! = Y2=X1*SIN(ang)+Y1*COS(ang)+Y0 =
! ===========================================================
IMPLICIT REAL(8)(A-H,O-Z)
REAL(8) N(6),U(3)
CHARACTER SIGN*1,NAME*10
DATA IFL3/4/
COMMON AVE_X,AVE_Y
SIGN='+'
PI=4.D0*DATAN(1.D0)
ANG=0.
X0=0.
Y0=0.
AVE_X=0.
AVE_Y=0.
REWIND IFL3
! === 求第一坐標系統之形心位置
DO I=1,NUM
READ(IFL3) NAME,X1,Y1,X2,Y2
AVE_X=AVE_X+X1
AVE_Y=AVE_Y+Y1
END DO
AVE_X=AVE_X/NUM
AVE_Y=AVE_Y/NUM
WRITE(*,900) AVE_Y,AVE_X
WRITE(2,900) AVE_Y,AVE_X
900 FORMAT(/10X,'旋轉中心坐標',2X,'AVE_N=',F13.3,5X,'AVE_E=',F13.3)
DO 40 ITER=1,20
! WRITE(*,1004) ITER
1004 FORMAT(/2X,'#### ITERATION ',I2,' #### ')
DO 10 I=1,6
N(I)=0.
10 CONTINUE
DO 20 I=1,3
U(I)=0.
20 CONTINUE
! === 法方程式 Normal Equations
REWIND IFL3
DO 30 I=1,NUM
READ(IFL3) NAME,X1,Y1,X2,Y2
A1=-(X1-AVE_X)*DSIN(ANG)-(Y1-AVE_Y)*DCOS(ANG)
A2= (X1-AVE_X)*DCOS(ANG)-(Y1-AVE_Y)*DSIN(ANG)
N(1)=N(1)+A1*A1+A2*A2
N(2)=N(2)+A1
N(3)=N(3)+1.0
N(4)=N(4)+A2
N(6)=N(6)+1.0
B1=(X2-AVE_X)-(X1-AVE_X)*DCOS(ANG)+(Y1-AVE_Y)*DSIN(ANG)-X0
B2=(Y2-AVE_Y)-(X1-AVE_X)*DSIN(ANG)-(Y1-AVE_Y)*DCOS(ANG)-Y0
U(1)=U(1)+A1*B1+A2*B2
U(2)=U(2)+B1
U(3)=U(3)+B2
30 CONTINUE
CALL SINV(N,3,1.D-9,IER)
! WRITE(*,1002) IER
1002 FORMAT(2X,'IER=',I2)
IF(IER.EQ.-1) THEN
WRITE(*,1200)
1200 FORMAT(//10X,'錯誤 : 法係數矩陣奇異解(singular) ! 程式終止 !')
STOP
END IF
! WRITE(*,1000) N(1),N(2),N(4),N(3),N(5),N(6)
1000 FORMAT(/5X,'INVERSE OF NORMAL EQ.'/&
8X,3(E14.8,2X)/24X,2(E14.8,2X)/40X,1(E14.8,2X))
COEF1=N(2)/DSQRT(N(1)*N(3))
COEF2=N(4)/DSQRT(N(1)*N(6))
COEF3=N(5)/DSQRT(N(3)*N(6))
! WRITE(*,1010) 1.,COEF1,COEF2,1.,COEF3,1.
1010 FORMAT(/5X,'法係數矩陣之相關係數 CORRELATION COEFFICIENT'/&
8X,3(E14.8,2X)/24X,2(E14.8,2X)/40X,1(E14.8,2X))
DANG=N(1)*U(1)+N(2)*U(2)+N(4)*U(3)
DX0 =N(2)*U(1)+N(3)*U(2)+N(5)*U(3)
DY0 =N(4)*U(1)+N(5)*U(2)+N(6)*U(3)
ANG=ANG+DANG
X0=X0+DX0
Y0=Y0+DY0
! IF(DABS(DANG).LE.1.D-6) GOTO 50
! IF((DABS(DANG).LE.(0.00001/206265.D0)).AND.(DABS(DX0).LE.0.0001)
! $ .AND.(DABS(DY0).LE.0.0001)) GOTO 50
40 CONTINUE
50 PVV=0.
REWIND IFL3
WRITE(*,1020)
WRITE(2,1020)
1020 FORMAT(/1x,'[控 制 點 轉 換 平 差 成 果]'/&
5X,'點 名',8X,'轉換前(N,E)坐標',9X,'轉換後(N,E)坐標',&
7X,'不符值(dN,dE)'/80(1H=))
DO 60 I=1,NUM
READ(IFL3) NAME,X1,Y1,X2,Y2
VX=(X1-AVE_X)*DCOS(ANG)-(Y1-AVE_Y)*DSIN(ANG)+X0-(X2-AVE_X)
VY=(X1-AVE_X)*DSIN(ANG)+(Y1-AVE_Y)*DCOS(ANG)+Y0-(Y2-AVE_Y)
WRITE(*,1030) I,NAME,Y1,X1,Y2,X2,VY,VX
WRITE(2,1030) I,NAME,Y1,X1,Y2,X2,VY,VX
1030 FORMAT(I3,2X,A10,2(2X,F11.3,1X,F10.3),2X,F7.3,1X,F7.3)
PVV=PVV+VX*VX+VY*VY
60 CONTINUE
IDF=2*NUM-3
IF(IDF.EQ.0) THEN
SIGMA=0.
ELSE
SIGMA=DSQRT(PVV/IDF)
SA=SIGMA*DSQRT(N(1))
SB=SIGMA*DSQRT(N(3))
SC=SIGMA*DSQRT(N(6))
TANG=ABS(ANG/SA)
TX0=ABS(X0/SB)
TY0=ABS(Y0/SC)
END IF
CALL DEGREE(ANG,IDEG,MIN,SEC)
IF(ANG.LT.0.) SIGN='-'
WRITE(*,1040) SIGN,ABS(IDEG),MIN,SEC,SA,TANG,Y0,SC,TY0,&
X0,SB,TX0,SIGMA,ITER-1,IDF
WRITE(2,1040) SIGN,ABS(IDEG),MIN,SEC,SA,TANG,Y0,SC,TY0,&
X0,SB,TX0,SIGMA,ITER-1,IDF
1040 FORMAT(/10X,'旋轉量 = ',A1,I3,'d',I2,'m',F8.5,'s',&
2X,'中誤差=',E10.4,2X,'t=',E8.2/&
10X,'N軸平移量= ',F17.4,2X,'中誤差=',E10.4,2X,'t=',E8.2/&
10X,'E軸平移量= ',F17.4,2X,'中誤差=',E10.4,2X,'t=',E8.2/&
10X,'單位權中誤差=',F8.5,9X,'漸近次數=',I3,7X,'自由度=',I5)
WRITE(2,1050) 1.,COEF1,COEF2,1.,COEF3,1.
1050 FORMAT(/10X,'法係數矩陣之相關係數Correlation Coefficient Matrix'/&
10x,'rotation',8X,3(E10.4,2X)/&
10x,'x0 ',20X,2(E10.4,2X)/&
10x,'y0 ',32X,1(E10.4,2X))
WRITE(8,*) AVE_Y
WRITE(8,*) AVE_X
WRITE(8,*) (ANG*180.d0/pi)
WRITE(8,*) X0
WRITE(8,*) Y0
RETURN
END
SUBROUTINE TR3(X,Y,X0,Y0,ANG,XX,YY)
IMPLICIT REAL(8) (A-H,O-Z)
COMMON AVE_X,AVE_Y
XX=(X-AVE_X)*DCOS(ANG)-(Y-AVE_Y)*DSIN(ANG)+X0+AVE_X
YY=(X-AVE_X)*DSIN(ANG)+(Y-AVE_Y)*DCOS(ANG)+Y0+AVE_Y
RETURN
END
SUBROUTINE DEGREE(RAD,IDEG,MIN,SEC)
! ================================================================
! = PROGRAM : TRANSFORM RADIAN TO DEGREE-MINUTE-SECOND =
! = modify PI 2023.5.12 =
! = 更正 IF(KEY.EQ.1) IDEG=-IDEG 位置 =
! ================================================================
IMPLICIT REAL(8) (A-H,O-Z)
! PI=0.31415926536D01
PI=4.D0*DATAN(1.D0)
KEY=0
DEG=RAD*180.D0/PI
IF(DEG.LT.0.) THEN
DEG=-DEG
KEY=1
END IF
IDEG=INT(DEG)
SMIN=(DEG-IDEG)*60.D0
MIN=INT(SMIN)
SEC=(SMIN-MIN)*60.D0
IF(KEY.EQ.1) IDEG=-IDEG
RETURN
END
! =========================================================
! = Modify 2000.12.30 固定輸出入檔案之檔名 =
! = 2001.5.24 修改計算模式 =
! = 90.05.29. 修改轉換參數之輸出有效位數(蕭欽霖之貢獻) =
! = 90.10.14. 輸出轉換參數值, 以供AutoCAD程式使用 =
! = 91.05.09. 修改不符值輸出順序錯誤 =
! = 2004.12.15 compiled by Lahey Fortran 95 for Windows XP =
! = 2014.01.28 將旋轉中心(0,0)移至測區中央 =
! = 適用於兩坐標系統均為二度分帶系統 =
! = 2023.02.19 修正度分秒副程式之bug =
! = 以GNU Fortran Compile =
! ========================================================
IMPLICIT REAL(8)(A-H,O-Z)
REAL(8) XF,YF,XI,YI
CHARACTER INFILE*60,OUTFILE*60,CORFILE*60,JOB*70,PARFILE*60
CHARACTER*10 NAME,PNAME
INTEGER idt(8)
CHARACTER*10 date,time,zone
COMMON AVEX1,AVEY1,AVEX2,AVEY2
! CALL GETDAT(IYEAR,IMONTH,IDAY)
! CALL GETTIM(IHOUR,IMIN,ISEC,I100H)
CALL date_and_time(date,time,zone,idt)
WRITE(*,1000)
1000 FORMAT(20X,'***************************************'/&
20X,'* HELMERT 三參數相似轉換 *'/&
20X,'* 形心對形心坐標轉換 *'/&
20X,'* compiler GNU Fortran 95 *'/&
20X,'* program by ChinaTom *'/&
20X,'* Version 3.0 2023.5. *'/&
20X,'***************************************')
WRITE(*,1002)
1002 FORMAT(/10X,'本程式適用於兩坐標系統均為二度分帶系統,且將旋轉中心'&
'(0,0)移至測區中央')
WRITE(*,1010)
1010 FORMAT(/10X,'輸入資料檔名 <<主檔名即可>>： '$)
READ(*,1020) INFILE
1020 FORMAT(A60)
NN=LEN_TRIM(INFILE)
ICHK=INDEX(INFILE,'.')
IF(ICHK.NE.0) THEN
WRITE(*,1030)
1030 FORMAT(/10X,'<< 輸入錯誤 : 檔名中含附加檔名 >> 程式中止!')
! PAUSE ' 程式中止! 按 ENTER 鍵離開'
STOP
END IF
OUTFILE=INFILE
CORFILE=INFILE
PARFILE=INFILE
INFILE(NN+1:NN+4)= '.DAT'
OUTFILE(NN+1:NN+4)='.DOC'
CORFILE(NN+1:NN+4)='.COD'
PARFILE(NN+1:NN+4)='.PA3'
OPEN(1,FILE=INFILE,STATUS='OLD',IOSTAT=ICHK)
! IF(ICHK.EQ.6416) THEN
IF(ICHK.GT.0) THEN
WRITE(*,1040) INFILE
1040 FORMAT(/10X,'<<錯誤>> : 找不到資料檔 ',A,' 程式中止!')
! PAUSE ' 程式中止! 按 ENTER 鍵離開'
STOP
END IF
OPEN(2,FILE=OUTFILE,STATUS='UNKNOWN')
OPEN(3,FILE=CORFILE,STATUS='UNKNOWN')
OPEN(4,FILE='AA.TMP',STATUS='UNKNOWN',FORM='UNFORMATTED')
OPEN(5,FILE='BB.TMP',STATUS='UNKNOWN',FORM='UNFORMATTED')
OPEN(8,FILE=PARFILE,STATUS='UNKNOWN')
WRITE(2,1000)
NUM=1
NUMRJ=0
READ(1,'(A70)') JOB
10 READ(1,1050,END=20) IC,NAME,YI,XI,YF,XF
1050 FORMAT(I5,A10,4F15.0)
! =============================================================
! = IC : >0 POINT NOT USED, =-99 END OF CONTROL POINTS =
! = NAME : POINT NAME =
! = XI,YI : (X,Y) OF ORIGINAL COORDS. SYSTEM =
! = XF,YF : (X,Y) OF DESTINATION COORDS. SYSTEM =
! =============================================================
IF(IC.GT.0) THEN
WRITE(5) NAME,XI,YI,XF,YF
NUMRJ=NUMRJ+1
GOTO 10
END IF
IF(IC.EQ.-99) GOTO 20
WRITE(4) NAME,XI,YI,XF,YF
NUM=NUM+1
GOTO 10
20 NUM=NUM-1
IF(NUM.LE.1) THEN
WRITE(*,1060)
1060 FORMAT(//10X,'錯誤：控制點點數少於 2點 !')
STOP
END IF
WRITE(2,1070) TRIM(JOB)
1070 FORMAT(/5X,'資料名稱 : ',A)
WRITE(2,1080) TRIM(INFILE)
1080 FORMAT( 5X,'檔案名稱 : ',A)
! WRITE(2,1090) IYEAR,IMONTH,IDAY,IHOUR,IMIN
!1090 FORMAT( 5X,'計算日期時間 : ',I4,'/',I2,'/',I2,2X,I2,':',I2)
WRITE(2,1090) IDT(1),IDT(2),IDT(3),IDT(5),IDT(6),IDT(7)
1090 FORMAT(5X,'計算日期 : ',I4,'/',I2,'/',I2,4X,'時間: ',I2,':',I2,':',I2)
IF(NUMRJ.NE.0) THEN
REWIND 5
WRITE(2,1100)
1100 FORMAT(/1X,'[未加入轉換之已知坐標點位]'/&
11X,'點 名',11X,'轉換前(N,E)坐標',12X,'轉換後(N,E)坐標'/&
3X,72(1H=))
DO 22 I=1,NUMRJ
READ(5) NAME,XI,YI,XF,YF
WRITE(2,1110) I,NAME,YI,XI,YF,XF
1110 FORMAT(4X,I3,3X,A10,2(3X,F11.3,2X,F11.3))
22 CONTINUE
END IF
CALL PARA3(X04,Y04,ANG4,NUM)
! IF(EOF(1)) GOTO 40
WRITE(*,1120)
WRITE(2,1120)
1120 FORMAT(//1X,'[未 知 點 轉 換 成 果]'/&
11X,'點 名',11X,'轉換前(N,E)坐標',12X,'轉換後(N,E)坐標'/&
8X,67(1H=))
30 READ(1,1130,END=40) IC,PNAME,Y,X
1130 FORMAT(I5,A10,2F15.0)
IF(IC.GT.0) GOTO 30
CALL TR3(X,Y,X04,Y04,ANG4,XX,YY)
WRITE(*,1140) PNAME,Y,X,YY,XX
WRITE(2,1140) PNAME,Y,X,YY,XX
1140 FORMAT(10X,A10,2(3X,F11.3,2X,F11.3))
WRITE(3,1150) PNAME,YY,XX
1150 FORMAT(A10,2F15.3)
GOTO 30
40 CONTINUE
CLOSE(4,STATUS='DELETE')
CLOSE(5,STATUS='DELETE')
CLOSE(8)
END
SUBROUTINE PARA3(X0,Y0,ANG,NUM)
! ===========================================================
! = PROGRAM : 3-PARAMETER TRANSFORMATION =
! = (X,Y)-SHIFT, ROTATION =
! =---------------------------------------------------------=
! = X2=X1*COS(ang)-Y1*SIN(ang)+X0 =
! = Y2=X1*SIN(ang)+Y1*COS(ang)+Y0 =
! ===========================================================
IMPLICIT REAL(8)(A-H,O-Z)
REAL(8) N(6),U(3)
CHARACTER SIGN*1,NAME*10
DATA IFL3/4/
COMMON AVE_X,AVE_Y
SIGN='+'
PI=4.D0*DATAN(1.D0)
ANG=0.
X0=0.
Y0=0.
AVE_X=0.
AVE_Y=0.
REWIND IFL3
! === 求第一坐標系統之形心位置
DO I=1,NUM
READ(IFL3) NAME,X1,Y1,X2,Y2
AVE_X=AVE_X+X1
AVE_Y=AVE_Y+Y1
END DO
AVE_X=AVE_X/NUM
AVE_Y=AVE_Y/NUM
WRITE(*,900) AVE_Y,AVE_X
WRITE(2,900) AVE_Y,AVE_X
900 FORMAT(/10X,'旋轉中心坐標',2X,'AVE_N=',F13.3,5X,'AVE_E=',F13.3)
DO 40 ITER=1,20
! WRITE(*,1004) ITER
1004 FORMAT(/2X,'#### ITERATION ',I2,' #### ')
DO 10 I=1,6
N(I)=0.
10 CONTINUE
DO 20 I=1,3
U(I)=0.
20 CONTINUE
! === 法方程式 Normal Equations
REWIND IFL3
DO 30 I=1,NUM
READ(IFL3) NAME,X1,Y1,X2,Y2
A1=-(X1-AVE_X)*DSIN(ANG)-(Y1-AVE_Y)*DCOS(ANG)
A2= (X1-AVE_X)*DCOS(ANG)-(Y1-AVE_Y)*DSIN(ANG)
N(1)=N(1)+A1*A1+A2*A2
N(2)=N(2)+A1
N(3)=N(3)+1.0
N(4)=N(4)+A2
N(6)=N(6)+1.0
B1=(X2-AVE_X)-(X1-AVE_X)*DCOS(ANG)+(Y1-AVE_Y)*DSIN(ANG)-X0
B2=(Y2-AVE_Y)-(X1-AVE_X)*DSIN(ANG)-(Y1-AVE_Y)*DCOS(ANG)-Y0
U(1)=U(1)+A1*B1+A2*B2
U(2)=U(2)+B1
U(3)=U(3)+B2
30 CONTINUE
CALL SINV(N,3,1.D-9,IER)
! WRITE(*,1002) IER
1002 FORMAT(2X,'IER=',I2)
IF(IER.EQ.-1) THEN
WRITE(*,1200)
1200 FORMAT(//10X,'錯誤 : 法係數矩陣奇異解(singular) ! 程式終止 !')
STOP
END IF
! WRITE(*,1000) N(1),N(2),N(4),N(3),N(5),N(6)
1000 FORMAT(/5X,'INVERSE OF NORMAL EQ.'/&
8X,3(E14.8,2X)/24X,2(E14.8,2X)/40X,1(E14.8,2X))
COEF1=N(2)/DSQRT(N(1)*N(3))
COEF2=N(4)/DSQRT(N(1)*N(6))
COEF3=N(5)/DSQRT(N(3)*N(6))
! WRITE(*,1010) 1.,COEF1,COEF2,1.,COEF3,1.
1010 FORMAT(/5X,'法係數矩陣之相關係數 CORRELATION COEFFICIENT'/&
8X,3(E14.8,2X)/24X,2(E14.8,2X)/40X,1(E14.8,2X))
DANG=N(1)*U(1)+N(2)*U(2)+N(4)*U(3)
DX0 =N(2)*U(1)+N(3)*U(2)+N(5)*U(3)
DY0 =N(4)*U(1)+N(5)*U(2)+N(6)*U(3)
ANG=ANG+DANG
X0=X0+DX0
Y0=Y0+DY0
! IF(DABS(DANG).LE.1.D-6) GOTO 50
! IF((DABS(DANG).LE.(0.00001/206265.D0)).AND.(DABS(DX0).LE.0.0001)
! $ .AND.(DABS(DY0).LE.0.0001)) GOTO 50
40 CONTINUE
50 PVV=0.
REWIND IFL3
WRITE(*,1020)
WRITE(2,1020)
1020 FORMAT(/1x,'[控 制 點 轉 換 平 差 成 果]'/&
5X,'點 名',8X,'轉換前(N,E)坐標',9X,'轉換後(N,E)坐標',&
7X,'不符值(dN,dE)'/80(1H=))
DO 60 I=1,NUM
READ(IFL3) NAME,X1,Y1,X2,Y2
VX=(X1-AVE_X)*DCOS(ANG)-(Y1-AVE_Y)*DSIN(ANG)+X0-(X2-AVE_X)
VY=(X1-AVE_X)*DSIN(ANG)+(Y1-AVE_Y)*DCOS(ANG)+Y0-(Y2-AVE_Y)
WRITE(*,1030) I,NAME,Y1,X1,Y2,X2,VY,VX
WRITE(2,1030) I,NAME,Y1,X1,Y2,X2,VY,VX
1030 FORMAT(I3,2X,A10,2(2X,F11.3,1X,F10.3),2X,F7.3,1X,F7.3)
PVV=PVV+VX*VX+VY*VY
60 CONTINUE
IDF=2*NUM-3
IF(IDF.EQ.0) THEN
SIGMA=0.
ELSE
SIGMA=DSQRT(PVV/IDF)
SA=SIGMA*DSQRT(N(1))
SB=SIGMA*DSQRT(N(3))
SC=SIGMA*DSQRT(N(6))
TANG=ABS(ANG/SA)
TX0=ABS(X0/SB)
TY0=ABS(Y0/SC)
END IF
CALL DEGREE(ANG,IDEG,MIN,SEC)
IF(ANG.LT.0.) SIGN='-'
WRITE(*,1040) SIGN,ABS(IDEG),MIN,SEC,SA,TANG,Y0,SC,TY0,&
X0,SB,TX0,SIGMA,ITER-1,IDF
WRITE(2,1040) SIGN,ABS(IDEG),MIN,SEC,SA,TANG,Y0,SC,TY0,&
X0,SB,TX0,SIGMA,ITER-1,IDF
1040 FORMAT(/10X,'旋轉量 = ',A1,I3,'d',I2,'m',F8.5,'s',&
2X,'中誤差=',E10.4,2X,'t=',E8.2/&
10X,'N軸平移量= ',F17.4,2X,'中誤差=',E10.4,2X,'t=',E8.2/&
10X,'E軸平移量= ',F17.4,2X,'中誤差=',E10.4,2X,'t=',E8.2/&
10X,'單位權中誤差=',F8.5,9X,'漸近次數=',I3,7X,'自由度=',I5)
WRITE(2,1050) 1.,COEF1,COEF2,1.,COEF3,1.
1050 FORMAT(/10X,'法係數矩陣之相關係數Correlation Coefficient Matrix'/&
10x,'rotation',8X,3(E10.4,2X)/&
10x,'x0 ',20X,2(E10.4,2X)/&
10x,'y0 ',32X,1(E10.4,2X))
WRITE(8,*) AVE_Y
WRITE(8,*) AVE_X
WRITE(8,*) (ANG*180.d0/pi)
WRITE(8,*) X0
WRITE(8,*) Y0
RETURN
END
SUBROUTINE TR3(X,Y,X0,Y0,ANG,XX,YY)
IMPLICIT REAL(8) (A-H,O-Z)
COMMON AVE_X,AVE_Y
XX=(X-AVE_X)*DCOS(ANG)-(Y-AVE_Y)*DSIN(ANG)+X0+AVE_X
YY=(X-AVE_X)*DSIN(ANG)+(Y-AVE_Y)*DCOS(ANG)+Y0+AVE_Y
RETURN
END
SUBROUTINE DEGREE(RAD,IDEG,MIN,SEC)
! ================================================================
! = PROGRAM : TRANSFORM RADIAN TO DEGREE-MINUTE-SECOND =
! = modify PI 2023.5.12 =
! = 更正 IF(KEY.EQ.1) IDEG=-IDEG 位置 =
! ================================================================
IMPLICIT REAL(8) (A-H,O-Z)
! PI=0.31415926536D01
PI=4.D0*DATAN(1.D0)
KEY=0
DEG=RAD*180.D0/PI
IF(DEG.LT.0.) THEN
DEG=-DEG
KEY=1
END IF
IDEG=INT(DEG)
SMIN=(DEG-IDEG)*60.D0
MIN=INT(SMIN)
SEC=(SMIN-MIN)*60.D0
IF(KEY.EQ.1) IDEG=-IDEG
RETURN
END
SUBROUTINE SINV(A,N,EPS,IER)
IMPLICIT REAL(8) (A-H,O-Z)
REAL(8) A(1),DIN,WORK,DPIV,DSUM
IF(N-1) 12,1,1
1 IER=0
KPIV=0
DO 11 K=1,N
KPIV=KPIV+K
IND=KPIV
LEND=K-1
TOL=DABS(EPS*(A(KPIV)))
DO 11 I=K,N
DSUM=0.D0
IF(LEND) 2,4,2
2 DO 3 L=1,LEND
LANF=KPIV-L
LIND=IND-L
3 DSUM=DSUM+(A(LANF)*A(LIND))
4 DSUM=A(IND)-DSUM
IF(I-K) 10,5,10
5 IF(DSUM-TOL) 6,6,9
6 IF(DSUM) 12,12,7
7 IF(IER) 8,8,9
8 IER=K-1
9 DPIV=DSQRT(DSUM)
A(KPIV)=DPIV
DPIV=1.0/DPIV
GO TO 11
10 A(IND)=DSUM*DPIV
11 IND=IND+I
GO TO 13
12 IER=-1
13 IF(IER) 22,14,14
14 IPIV=N*(N+1)/2
IND=IPIV
DO 19 I=1,N
DIN=1.0/A(IPIV)
A(IPIV)=DIN
MIN=N
KEND=I-1
LANF=N-KEND
IF(KEND) 18,18,15
15 J=IND
DO 17 K=1,KEND
WORK=0.0
MIN=MIN-1
LHOR=IPIV
LVER=J
DO 16 L=LANF,MIN
LVER=LVER+1
LHOR=LHOR+L
16 WORK=WORK+(A(LVER)*A(LHOR))
A(J)=-WORK*DIN
17 J=J-MIN
18 IPIV=IPIV-MIN
19 IND=IND-1
DO 21 I=1,N
IPIV=IPIV+I
J=IPIV
DO 21 K=I,N
WORK=0.0
LHOR=J
DO 20 L=K,N
LVER=LHOR+K-I
WORK=WORK+(A(LHOR)*A(LVER))
20 LHOR=LHOR+L
A(J)=WORK
21 J=J+K
22 RETURN
END
IMPLICIT REAL(8) (A-H,O-Z)
REAL(8) A(1),DIN,WORK,DPIV,DSUM
IF(N-1) 12,1,1
1 IER=0
KPIV=0
DO 11 K=1,N
KPIV=KPIV+K
IND=KPIV
LEND=K-1
TOL=DABS(EPS*(A(KPIV)))
DO 11 I=K,N
DSUM=0.D0
IF(LEND) 2,4,2
2 DO 3 L=1,LEND
LANF=KPIV-L
LIND=IND-L
3 DSUM=DSUM+(A(LANF)*A(LIND))
4 DSUM=A(IND)-DSUM
IF(I-K) 10,5,10
5 IF(DSUM-TOL) 6,6,9
6 IF(DSUM) 12,12,7
7 IF(IER) 8,8,9
8 IER=K-1
9 DPIV=DSQRT(DSUM)
A(KPIV)=DPIV
DPIV=1.0/DPIV
GO TO 11
10 A(IND)=DSUM*DPIV
11 IND=IND+I
GO TO 13
12 IER=-1
13 IF(IER) 22,14,14
14 IPIV=N*(N+1)/2
IND=IPIV
DO 19 I=1,N
DIN=1.0/A(IPIV)
A(IPIV)=DIN
MIN=N
KEND=I-1
LANF=N-KEND
IF(KEND) 18,18,15
15 J=IND
DO 17 K=1,KEND
WORK=0.0
MIN=MIN-1
LHOR=IPIV
LVER=J
DO 16 L=LANF,MIN
LVER=LVER+1
LHOR=LHOR+L
16 WORK=WORK+(A(LVER)*A(LHOR))
A(J)=-WORK*DIN
17 J=J-MIN
18 IPIV=IPIV-MIN
19 IND=IND-1
DO 21 I=1,N
IPIV=IPIV+I
J=IPIV
DO 21 K=I,N
WORK=0.0
LHOR=J
DO 20 L=K,N
LVER=LHOR+K-I
WORK=WORK+(A(LHOR)*A(LVER))
20 LHOR=LHOR+L
A(J)=WORK
21 J=J+K
22 RETURN
END
留言
張貼留言