program stalk_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 funcs integer b parameter(b=2) integer ma,mfit,ncvm,ndata,lista(b),ii integer nmbr_stalks real sig(5),z(5),d(5) real a(b),chisq,covar(b,b),lngth real inch_to_cm real a1_ave,a2_ave open(unit=21,file='stalk_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/ print*,'ok' a1_ave=0.0 a2_ave=0.0 read(21,*) nmbr_stalks print*,'nmbr_stalks',nmbr_stalks do ii=1,nmbr_stalks lista(1)=1 lista(2)=2 sig(1)=1.0 sig(2)=1.0 sig(3)=1.0 sig(4)=1.0 sig(5)=1.0 inch_to_cm=2.54 ma=b ncvm=b mfit=b ndata=5 ccccccccccccccccccccccccccccccccccccccccccccccc c Read measured (polar) data ccccccccccccccccccccccccccccccccccccccccccccccc read(21,*) z(1),d(1) read(21,*) z(2),d(2) read(21,*) z(3),d(3) read(21,*) z(4),d(4) read(21,*) z(5),d(5) read(21,*) lngth lngth=lngth*inch_to_cm ccccccccccccccccccccccccccccccccccccccccccccccc c Convert to cartesian data and write c to "geom.dat" ccccccccccccccccccccccccccccccccccccccccccccccc z(1)=z(1)*inch_to_cm/lngth z(2)=z(2)*inch_to_cm/lngth z(3)=z(3)*inch_to_cm/lngth z(4)=z(4)*inch_to_cm/lngth z(5)=z(5)*inch_to_cm/lngth d(1)=d(1)*inch_to_cm d(2)=d(2)*inch_to_cm d(3)=d(3)*inch_to_cm d(4)=d(4)*inch_to_cm d(5)=d(5)*inch_to_cm ccccccccccccccccccccccccccccccccccccccccccccccc c Determine basis function coefficients cccccccccccccccccccccccccccccccccccccccccccccc print*,'ok',ii call lfit(z,d,sig,ndata,a,ma,lista,mfit, & covar,ncvm,chisq,funcs) cccccccccccccccccccccccccccccccccccccccccccccccccc c Write coefficients and chisq to "geom.out" cccccccccccccccccccccccccccccccccccccccccccccccccc print*, a(1),a(2),chisq,lngth a1_ave=a1_ave+a(1)/nmbr_stalks a2_ave=a2_ave+a(2)/nmbr_stalks enddo close(21) print*, a1_ave,a2_ave/a1_ave 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