program geom ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c version 1.0 8/19/93 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c This program is written to provide a simple polynomial estimate c of the shape of a wheat leaf blade. Four points along a leaf c blade where measured in the field (polar coordinates) and these c points are contained in "geom.in". A set of basis functions are c then selected, as specified in subroutine "funcs". Variable "b" c is set to the number of basis functions. Subroutine lfit (et al) c is then called to determine the real coefficients of the basis c functions such that the sum of these weighted functions provides c the minimum mean-squared error between the estimate and the c measured data. The coefficients are written to file "geom.out", c along with "chisq", a measurement of the estimate error. c The measured data, converted to cartesian coordinates, are c written to "geom.dat". Both "geom.out" and "geom.dat" are c input files to "test.ftn", which plots the results of this c basis function estimate along with the measured data. c c NOTE: PROGRAM MUST BE BOUND WITH SUB "FUNCS.BIN" ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc external funncs integer b parameter(b=2) integer ma,mfit,ncvm,ndata,lista(b),i real sig(4),m(4),p(4),x(4),y(4),pi,deg real a(b),chisq,covar(b,b),lngth open(unit=20,file='geom.out') open(unit=21,file='geom.dat') ccccccccccccccccccccccccccccccccccccccccccccccc c Initialize data ccccccccccccccccccccccccccccccccccccccccccccccc c ATTENTION! the data list of the next line must proceed c to the value of b. (example: b=3 then -> data lista /1,2,3/ data lista /1,2/ data sig /4 * 1.0/ pi=4.0*atan(1.0) deg=pi/180. ma=b ncvm=b mfit=b ndata=4 ccccccccccccccccccccccccccccccccccccccccccccccc c Read measured (polar) data ccccccccccccccccccccccccccccccccccccccccccccccc read *, m(1), p(1) read *, m(2), p(2) read *, m(3), p(3) read *, m(4), p(4) read *, lngth ccccccccccccccccccccccccccccccccccccccccccccccc c Convert to cartesian data and write c to "geom.dat" ccccccccccccccccccccccccccccccccccccccccccccccc do 10 i=1,4 x(i)=m(i)*sin(p(i)*deg) y(i)=m(i)*cos(p(i)*deg) write(21,*) x(i),y(i) 10 continue ccccccccccccccccccccccccccccccccccccccccccccccc c Determine basis function coefficients cccccccccccccccccccccccccccccccccccccccccccccc call lfit(x,y,sig,ndata,a,ma,lista,mfit, & covar,ncvm,chisq,funncs) cccccccccccccccccccccccccccccccccccccccccccccccccc c Write coefficients and chisq to "geom.out" cccccccccccccccccccccccccccccccccccccccccccccccccc write(20,*) (a(j), j= 1,b),chisq,lngth close(20) close(21) stop end SUBROUTINE LFIT(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,COVAR,NCVM,CHISQ, *FUNCS) PARAMETER (MMAX=50) DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA), * COVAR(NCVM,NCVM),BETA(MMAX),AFUNC(MMAX) KK=MFIT+1 DO 12 J=1,MA IHIT=0 DO 11 K=1,MFIT IF (LISTA(K).EQ.J) IHIT=IHIT+1 11 CONTINUE IF (IHIT.EQ.0) THEN LISTA(KK)=J KK=KK+1 ELSE IF (IHIT.GT.1) THEN PAUSE 'Improper set in LISTA' ENDIF 12 CONTINUE IF (KK.NE.(MA+1)) PAUSE 'Improper set in LISTA' DO 14 J=1,MFIT DO 13 K=1,MFIT COVAR(J,K)=0. 13 CONTINUE BETA(J)=0. 14 CONTINUE DO 18 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) YM=Y(I) IF(MFIT.LT.MA) THEN DO 15 J=MFIT+1,MA YM=YM-A(LISTA(J))*AFUNC(LISTA(J)) 15 CONTINUE ENDIF SIG2I=1./SIG(I)**2 DO 17 J=1,MFIT WT=AFUNC(LISTA(J))*SIG2I DO 16 K=1,J COVAR(J,K)=COVAR(J,K)+WT*AFUNC(LISTA(K)) 16 CONTINUE BETA(J)=BETA(J)+YM*WT 17 CONTINUE 18 CONTINUE IF (MFIT.GT.1) THEN DO 21 J=2,MFIT DO 19 K=1,J-1 COVAR(K,J)=COVAR(J,K) 19 CONTINUE 21 CONTINUE ENDIF CALL GAUSSJ(COVAR,MFIT,NCVM,BETA,1,1) DO 22 J=1,MFIT A(LISTA(J))=BETA(J) 22 CONTINUE CHISQ=0. DO 24 I=1,NDATA CALL FUNCS(X(I),AFUNC,MA) SUM=0. DO 23 J=1,MA SUM=SUM+A(J)*AFUNC(J) 23 CONTINUE CHISQ=CHISQ+((Y(I)-SUM)/SIG(I))**2 24 CONTINUE CALL COVSRT(COVAR,NCVM,MA,LISTA,MFIT) RETURN END SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT) DIMENSION COVAR(NCVM,NCVM),LISTA(MFIT) DO 12 J=1,MA-1 DO 11 I=J+1,MA COVAR(I,J)=0. 11 CONTINUE 12 CONTINUE DO 14 I=1,MFIT-1 DO 13 J=I+1,MFIT IF(LISTA(J).GT.LISTA(I)) THEN COVAR(LISTA(J),LISTA(I))=COVAR(I,J) ELSE COVAR(LISTA(I),LISTA(J))=COVAR(I,J) ENDIF 13 CONTINUE 14 CONTINUE SWAP=COVAR(1,1) DO 15 J=1,MA COVAR(1,J)=COVAR(J,J) COVAR(J,J)=0. 15 CONTINUE COVAR(LISTA(1),LISTA(1))=SWAP DO 16 J=2,MFIT COVAR(LISTA(J),LISTA(J))=COVAR(1,J) 16 CONTINUE DO 18 J=2,MA DO 17 I=1,J-1 COVAR(I,J)=COVAR(J,I) 17 CONTINUE 18 CONTINUE RETURN END SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) PARAMETER (NMAX=50) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) DO 11 J=1,N IPIV(J)=0 11 CONTINUE DO 22 I=1,N BIG=0. DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'Singular matrix' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.' PIVINV=1./A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0. DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE RETURN END subroutine funncs(xx,aafunc,mma) c This subroutimne specifies the mma basis c functions used for programs geom.ftn and c test.ftn. integer mma real xx,aafunc(mma) c Specify basis functions aafunc(1)=xx aafunc(2)=xx**2 c aafunc(3)=xx**3 return end