real*4 mat(10000,10),YV(10000),p(25),ibad(100),iibad(100) real*8 mtm(10,10),ytm(10) INTEGER*2 H,K,L,HMAX,HOFF,HH(50000),KK(50000),LL(50000) INTEGER*2 I1(10,10,10),I2(200),ILINE(200) INTEGER*4 N,II,COUNT,UNIQUE,INDEX(0:50,0:50,0:70,1:5),INDEX2(50000,5) real*4 INT(50000),SIG(50000),XF(50000),YF(50000),FRAME(50000),PHI(100) REAL*4 CHI(50000),PSI(50000) C parameters: 20000 UNIQUE REFL, 100 FRAMES, COEF FOR UP TO 25 TERMS: C COEF(I,J,0) IS FOR COS(I*PHI+J*THETA), COEF(I,J,1) FOR SIN() real*4 IMEAN(20000),G(100),COEF(0:5,0:5,0:1) CHARACTER*40 FNAME,ASTRING CHARACTER*1 BLINE(400) EQUIVALENCE (ILINE,BLINE) R32K=128*256 pi=4.*atan(1.) write(6,*)' pi=',pi TYPE *,'READ DENZO OUTPUT FILE INTO ARRAY BY H,K,L;' TYPE *,' what was camera length?' read(5,*) dist TYPE *,' what are center coordinate in x and y (pixels)?' read(5,*) xcent,ycent TYPE *,' what are pixel dimensions in x and y dimension?' read(5,*) xpix,ypix TYPE *,' what ratio of I-Imean/I for rejection?' read(5,*) rratio 20 TYPE *,'CUMULATIVE SPOTS:',JC TYPE *,'Read (D)enzo output or (S)CALE LOCAL, or (W)rite sum?' READ(5,40)ASTRING IF (ASTRING.EQ.'D') GOTO 30 IF (ASTRING.EQ.'S') GOTO 100 IF (ASTRING.EQ.'W') GOTO 100 30 TYPE *,'CUMULATIVE SPOTS:',COUNT,'. ENTER FILE NAME (Q TO QUIT):' 35 read(5,40)FNAME IF (FNAME.EQ.'Q') GOTO 20 40 FORMAT (A20) OPEN(UNIT=2,FILE=FNAME,STATUS='OLD',ERR=500) write(6,*)'what is phi angle?' read(5,*)phif NFRAMES=NFRAMES+1 PHI(NFRAMES)=PHIF*pi/180. write(6,*)'phi=',phif,phi(nframes),nframes DO 45 I=1,5 45 READ (2,50)ASTRING C TYPE*, ASTRING 50 FORMAT (A40) N=0 55 READ (2,60)X 60 FORMAT (F7.4) 70 READ(2,101,END=500)H,K,L,J,X1,X2,X,S,X,X,Y if (H.EQ.999) GOTO 500 C TYPE *,'HKL=',H,K,L C 24 -8 71 1 66.6 66.6 4.64 0.0 0.971 35.4 623.0 0.034 0.0 C 28 -8 66 1 73.1 73.1 4.64 0.0 0.971 34.7 567.0 0.041 0.0 101 FORMAT (I4,2I4,I2,F8.0,F8.1,F7.2,F6.1,F6.3,F7.1,F7.1,F6.3,F8.1) IF (J.EQ.1) GOTO 70 if (x2.le.0) goto 70 COUNT=COUNT+1 IF (X1.GT.32000)TYPE*,H,K,L,X1 74 HH(COUNT)=H KK(COUNT)=K LL(COUNT)=L INT(COUNT)=X2 SIG(COUNT)=S XF(COUNT)=(X-xcent)*xpix YF(COUNT)=(Y-ycent)*ypix FRAME(COUNT)=NFRAMES GOTO 70 500 CLOSE (UNIT=2) C501 TYPE *,'NUMBER OF SPOTS REPORTED:',COUNT GOTO 30 C------------------------------------ 100 DO 120 II=1,COUNT H=IABS(HH(II)) K=IABS(KK(II)) L=IABS(LL(II)) C WRITE(6,*) II,H,K,L I=INDEX(H,K,L,1)+1 c write(6,*)' II, I=',ii,i IF (I.GT.5) THEN WRITE(6,*)' SKIPPING 5TH OBS FOR HKL=',H,K,L WRITE(6,*) HH(II),KK(II),LL(II),INT(II),SIG(II), . XF(II),YF(II),frame(II) do 110 j=2,5 jj=index(h,k,l,j) 110 write(6,*)hh(jj),kk(jj),ll(jj),frame(jj) GOTO 120 ENDIF IF (I.EQ.1) UNIQUE=UNIQUE+1 INDEX(H,K,L,1)=I INDEX(H,K,L,I+1)=II 120 CONTINUE WRITE(6,*) ' NUMBER UNIQUE REFL:',UNIQUE C------------------------------------- C COULD HAVE BEEN SETTING THIS UP WHEN MAKING INDEX: STORE THERE ONLY C REFERENCE TO INDEX2, PUT THE DATA IN INDEX2. BUT DON'T KNOW WHICH C WILL HAVE MULTIPLE MEASUREMENTS AT THAT TIME, SO HAVE TO SQUEZE C INDEX2 TO ELIMINATE SINGLE OBSERVATIONS. N=0 DO 200 H=1,50 DO 200 K=1,50 DO 200 L=1,50 IF (INDEX(H,K,L,1).GT.1) THEN N=N+1 DO 190 J=1,5 190 INDEX2(N,J)=INDEX(H,K,L,J) ENDIF 200 CONTINUE WRITE(6,*)'NUMBER OF REFL W MULTIPLE OBS:',N NMULT=N C----------------------------------------- C INITIALIZE ALL SCALE FACTORS TO 1.0: DO 210 I=1,NFRAMES 210 G(I)=1.0 C INITIALIZE COEF TO 0: DO 215 I=0,1 DO 215 J=0,1 DO 215 K=0,1 215 COEF(I,J,K)=0 coef(0,0,0)=1 C---------------------------------------- C CALCULATE CHI,PSI VALUES FOR EACH REFLECTION W MULT OBS DO 250 I=1,NMULT M=INDEX2(I,1) DO 250 J=1,M K=INDEX2(I,J+1) PHIF=PHI(FRAME(K)) CHI(K)=PHIF+ATAN(YF(K)/DIST) PSI(K)= ATAN(XF(K)/DIST) c write(6,*)'k,phif,chi(k),psi(k)=',k,phif,chi(k),psi(k) 250 continue C------------------------------------- C CALCULATE MEAN VALUES FOR EACH UNIQUE REFL: 280 continue c keep first frame scale factor = 1. gg=g(1) do 282 i=1,nframes 282 g(i)=g(i)/gg DO 300 II=1,NMULT M=INDEX2(II,1) SUMIG=0. SUMGG=0. DO 290 I=1,M J=INDEX2(II,1+I) GG=G(FRAME(J)) CH=CHI(J) PS=PSI(J) AB=0 DO 285 JJ=0,1 DO 285 MM=0,1 AB=AB+COEF(JJ,MM,0)*(COS(JJ*CH+MM*PS)) 285 AB=AB+COEF(JJ,MM,1)*(SIN(JJ*CH+MM*PS)) c WRITE(6,*)INT(J) SUMIG=SUMIG+ INT(J)*GG*AB 290 SUMGG=SUMGG+ GG*AB*GG*AB IMEAN(II)=SUMIG/SUMGG c WRITE(6,*)II,IMEAN(II) 300 CONTINUE C--------------------------------------- C CALCULATE DEVIATION FROM MEAN CUMN=0 CUMERR=0. nbad=0 DO 320 II=1,NMULT M=INDEX2(II,1) SUMIG=0. SUMGG=0. RMSERR=0. N=0 DO 310 I=1,M J=INDEX2(II,1+I) GG=G(FRAME(J)) CH=CHI(J) PS=PSI(J) AB=0 DO 305 JJ=0,1 DO 305 MM=0,1 AB=AB+COEF(JJ,MM,0)*(COS(JJ*CH+MM*PS)) 305 AB=AB+COEF(JJ,MM,1)*(SIN(JJ*CH+MM*PS)) c c WRITE(6,*)INT(J),INT(J)/(GG*AB),IMEAN(II) c ERR=(INT(J)/(GG*AB) - IMEAN(II)) ERR=(INT(J) - (GG*AB)*IMEAN(II)) if (abs(err/int(j)).gt.rratio) then nbad=nbad+1 ibad(nbad)=i iibad(nbad)=ii endif RMSERR=RMSERR+ERR*ERR N=N+1 SUMIG=SUMIG+ABS(ERR) 310 SUMGG=SUMGG+ABS(IMEAN(II)) c WRITE(6,*)II,'SUMOFI-IMEAN, SUMOFIMEAN, RATIO:', c . SUMIG,SUMGG,SUMIG/SUMGG c WRITE(6,*)'RMSERR',SQRT(RMSERR/N) CUMERR=CUMERR+RMSERR CUMN=CUMN+N 320 CONTINUE WRITE(6,*)'OVERALL RMSERR:',SQRT(CUMERR/CUMN) if (nbad.gt.0) then do 330 jj=nbad,1,-1 i=ibad(jj) ii=iibad(jj) m=index2(ii,1) write(6,*) ii,' bad:',(int(index2(ii,1+ij)),ij=1,m) do 325 j=i,m-1 325 index2(ii,j+1)=index2(ii,j+2) !jth index is in j+1 position! c so copying to jth from j+1th m=m-1 index2(ii,1)=m if(m.lt.1) then write(6,*) ii,' copying down rest of heap' do 328 k=ii,nmult-1 m=index2(k+1,1) do 328 j =1,m+1 328 index2(k,j)=index2(k+1,j) nmult=nmult-1 endif 330 continue do 340 ii=nmult,1,-1 if (index2(ii,1).lt.2) then write(6,*) ii,' copying down rest of heap' do 338 k=ii,nmult-1 m=index2(k+1,1) do 338 j =1, m+1 338 index2(k,j)=index2(k+1,j) nmult=nmult-1 endif 340 continue write(6,*)' nbad=', nbad write(6,*)' recalculating error after eliminating outliers.' goto 280 endif C------------------------------------- C CALCULATE FRAME SCALE FACTORS FOR EACH FRAME DO 400 II=1,NFRAMES SUMIG=0. SUMII=0. N=0 DO 390 I=1,NMULT M=INDEX2(I,1) DO 390 J=1,M K=INDEX2(I,J+1) IF (FRAME(K).EQ.II) THEN CH=CHI(K) PS=PSI(K) AB=0 DO 385 JJ=0,1 DO 385 MM=0,1 AB=AB+COEF(JJ,MM,0)*(COS(JJ*CH+MM*PS)) 385 AB=AB+COEF(JJ,MM,1)*(SIN(JJ*CH+MM*PS)) SUMII=SUMII+IMEAN(I)*AB*IMEAN(I)*AB SUMIG=SUMIG+INT(K)*IMEAN(I)*AB N=N+1 ENDIF 390 CONTINUE G(II)=SUMIG/SUMII 400 WRITE(6,*) 'FRAME ',II,' G=',G(II),' N=',N C------------------------------------ C CALCULATE COEFFICIENTS FOR ABSORPTION TERMS JJ=0 DO 552 II=1,NMULT M=INDEX2(II,1) DO 552 I=1,M JJ = JJ + 1 K=INDEX2(II,I+1) ch=chi(k) ps=psi(k) GI=IMEAN(II)*G(FRAME(K)) YV(JJ)=INT(K)-GI C MM=-1 DO 550 J=0,1 DO 550 L=0,1 MM=L+2*(J) C MM=MM+1 if (mm.eq.0) goto 550 c write(6,*)mm,mm+3,j,l MAT(jj,MM)=GI*(COS(J*CH+L*PS)) MAT(JJ,MM+3)=GI*(SIN(J*CH+L*PS)) 550 continue c write(6,553)yv(jj),(mat(jj,l1),l1=1,6) 552 continue 551 format(9e9.2) 553 format(9f9.0) M=JJ N=6 C LEAST SQUARES FITTING: c M is the number of observations, number of rows in mat() c N is the number of variables, number of columns in mat() c form i,j th element of mtm by dotting ith into j'th column of mat. c form ytm by dotting y unto each column of mat c multiply ytm by inverse of mtm, gives best coefficients. c form i,j th element of mtm by dotting ith into j'th column of mat. do 600 i=1,n do 600 j=1,i xdp=0 do 590 k=1,m 590 xdp=xdp+mat(k,i)*mat(k,j) mtm(i,j)=xdp 600 mtm(j,i)=xdp c form ytm by dotting y unto each column of mat do 700 i=1,n xdp=0 do 690 k=1,m 690 xdp=xdp+mat(k,i)*yv(k) 700 ytm(i)=xdp write(6,*)' matrix' do 730 ii=1,n 730 write(6,551) (mtm(ii,l),l=1,n) c multiply ytm by inverse of mtm, gives best coefficients. call matinv(mtm,n,x) write(6,*)' inverse matrix' do 740 ii=1,n 740 write(6,551) (mtm(ii,l),l=1,n) do 800 i=1,n x=0 do 790 j=1,n 790 x=x + mtm(i,j)*ytm(j) 800 p(i)=x C MM=0 DO 820 J=0,1 DO 820 L=0,1 MM=L+2*(J ) C MM=MM+1 if (mm.le.0) goto 820 write(6,*)j,l,p(mm),p(mm+3) COEF(J,L,0)=P(MM) COEF(J,L,1)=P(MM+3) 820 continue COEF(0,0,1)=0 write(6,*)'do another cycle?(y/n)' read(5,40)astring if (astring(:1).eq.'y') GOTO 280 end SUBROUTINE MATINV (ARRAY, NORDER, DET) DOUBLE PRECISION array, AMAX, SAVE 2 DIMENSION array(10,10), IK (10), JK (10) 10 DET = 1 11 DO 100 K=1, NORDER C C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX C AMAX=0 21 DO 30 I=K, NORDER DO 30 J=K, NORDER 23 IF (DABS (AMAX) - DABS (ARRAY(I,J) ) ) 24, 24, 30 24 AMAX=ARRAY(I,J) IK(K) =I JK(K) =J 30 CONTINUE C C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) C 31 IF (AMAX) 41, 32, 41 32 DET=0. GO TO 140 41 I=IK (K) IF (I-K) 21, 51, 43 43 DO 50 J=1, NORDER SAVE= ARRAY(K,J) ARRAY(K,J)= ARRAY(I,J) 50 ARRAY(I,J)= -SAVE 51 J=JK (K) IF (J-K) 21, 61, 53 53 DO 60 I=1, NORDER SAVE= ARRAY(I,K) ARRAY(I,K)= ARRAY(I,J) 60 ARRAY(I,J)= -SAVE C C ACCUMULATE ELEMENTS OF INVERSE MATRIX C 61 DO 70 I=1, NORDER IF(I-K) 63, 70, 63 63 ARRAY(I,K) = -ARRAY(I,K) / AMAX 70 CONTINUE 71 DO 80 I=1, NORDER DO 80 J=1, NORDER IF (I-K) 74, 80, 74 74 IF (J-K) 75, 80, 75 75 ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K) *ARRAY(K,J) 80 CONTINUE 81 DO 90 J=1, NORDER IF (J-K) 83, 90, 83 83 ARRAY(K,J) = ARRAY(K,J) / AMAX 90 CONTINUE ARRAY(K,K) = 1. / AMAX 100 continue c DET=DET *AMAX C C RESTORE ORDERING OF MATRIX C 101 DO 130 L=1, NORDER K= NORDER - L +1 J= IK (K) IF (J-K) 111, 111, 105 105 DO 110 I=1, NORDER SAVE= ARRAY(I,K) ARRAY(I,K)= -ARRAY(I,J) 110 ARRAY(I,J)= SAVE 111 I= JK(K) IF (I-K) 130, 130, 113 113 DO 120 J=1, NORDER SAVE= ARRAY(K,J) ARRAY(K,J) = -ARRAY(I,J) 120 ARRAY(I,J)= SAVE 130 CONTINUE 140 RETURN END