C C C MARXER.FOR C* ********************************************** C CHARACTER*80 TITLE CHARACTER*50 INFILE,OUTFILE,FILENAME CHARACTER*40 FNAME,ASTRING CHARACTER*1 BLINE(4000),OUTSTRING(4000),HSTRING(16),BOXFLAG INTEGER*2 ILINE(2000),H,K,L,N EQUIVALENCE (BLINE,ILINE) LOGICAL EFLG CHARACTER*12 AQDATE, AQTIME CHARACTER DAT*12,TIM*12 EQUIVALENCE (AQDATE,BLINE(101)),(AQTIME,BLINE(113)) INTEGER*4 xypix(2),lrecl, ICCOUNTS(2) real*4 xycenter(2) EQUIVALENCE (xypix,bline(1)),(lrecl, bline(9)) EQUIVALENCE (xycenter,bline(69)),(ICCOUNTS,BLINE(25)) INTEGER*4 EXPTIME EQUIVALENCE (EXPTIME,BLINE(33)) REAL*4 LAMBDA,DISTANCE,PHI_START,PHI_END EQUIVALENCE (LAMBDA,BLINE(77)),(DISTANCE,BLINE(81)),(PHI_START, - BLINE(85)),(PHI_END,BLINE(89)) CHARACTER*12 IPSYSTEM EQUIVALENCE (IPSYSTEM,BLINE(125)) R64K=256*256 R32K=128*256 DATA HSTRING/'0','1','2','3','4','5','6','7','8','9','A','B','C', - 'D','E','F'/ WRITE(6,1000) 1000 FORMAT(//' MAR RESEARCH TO POSTSCRIPT IMAGE FILE CONVERSION . PROGRAM Vers 950107 big plate' . /' To print image: LPR MARXER8.PS') 5 WRITE(6,1100) 1100 FORMAT(/'$Input filename (CR to quit): ') READ(5,1200) INFILE 1200 FORMAT(A) IF (INFILE(:5).EQ.' ') STOP INQUIRE(FILE=INFILE, NAME=FILENAME,EXIST=EFLG) IF (.NOT.EFLG) GOTO 5 OPEN (UNIT=2,FILE=INFILE,readonly,STATUS='OLD') TYPE *,'Enter source for coordinates of boxes around spots:' TYPE *,'CR-no boxes. D-denzo outfile. S-spots file from IMSTILLS' TYPE *,'O-oscgen generate file' READ(5,'(A1)')BOXFLAG IF (BOXFLAG.EQ.' ') GOTO 1400 IF (BOXFLAG.EQ.'d') BOXFLAG='D' IF (BOXFLAG.EQ.'s') BOXFLAG='S' IF (BOXFLAG.EQ.'o') BOXFLAG='O' IF (BOXFLAG.EQ.'S') STOP 'USE MX5 FOR NOW!' IF (BOXFLAG.EQ.'O') STOP 'USE MX6 FOR NOW!' 330 TYPE *,'ENTER FILENAME FOR BOX COORD FILE:' 340 FORMAT (A40) 335 read(5,340)FNAME 1400 write(6,*)'Enter window (pixels; xmin,xmax,ymin,ymax;', . ' 0 0 0 0 for default):' read(5,*) i1,i2,j1,j2 WRITE(6,*)'Enter density limits (min,max) or 0,0 to autoscale' WRITE(6,*)'Density autoscaling requires two passes through the file.' READ(5,*)X0,XM C WRITE(6,1600) C1600 FORMAT('$Output filename: ') C READ(5,1200) OUTFILE c OPEN(UNIT=3,FILE=OUTFILE,STATUS='NEW',RECL=1200) OPEN (UNIT=3,FILE='MARXER8.PS',STATUS='UNKNOWN',RECL=1200, & CARRIAGECONTROL='LIST') WRITE(3,'(A)')'%!PS-Adobe-3.0 ' WRITE(3,'(A)')'%%Creator: mx10' WRITE(3,'(A)')'%%Pages: 1 ' WRITE(3,'(A)')'%%BoundingBox:0 0 602 792 ' WRITE(3,'(A)')'%%EndComments ' WRITE(3,'(A)')'%%BeginProlog ' WRITE(3,'(A)')'%%EndProlog ' WRITE(3,'(A)')'%%Page: 1 1 ' WRITE(3,'(A)')'%%PageBoundingBox: 0 0 602 792 ' WRITE(3,'(A)')'initgraphics' c WRITE(3,'(A)')'%!PS-Adobe-3.0 EPSF-3.0 '//char(4) c WRITE(3,'(A)')'%%BoundingBox:0 0 602 792' write(3,'(A)')'/Times-Roman findfont 9 scalefont setfont' write(3,'(A)') & '/cshow {dup stringwidth pop 2 div neg 0 rmoveto show} def' CALL DATE(DAT) CALL TIME(TIM) write(3,'(A)') '50 772 moveto'// & '(Image File: '//FILENAME//' Printed:'//DAT//TIM//') show' C FIND MIN, MAX, MEAN DENSITY DMEAN=0 DMAX=-1E32 DMIN=1E32 C read HEADER into bline, encode into titles: READ(2,'(124A1)') (BLINE(II),II=1,124) 155 FORMAT (4096A1) call byteswapi4(xypix(1)) call byteswapi4(xypix(2)) call byteswapi4(lrecl) call byteswapi4(iccounts(1)) call byteswapi4(iccounts(2)) call byteswapr4(xycenter(1)) call byteswapr4(xycenter(2)) c npix=xypix(1) nbytes=lrecl if (nbytes.eq.1200) nbytes = 2400 NX=xypix(1) NY=xypix(2) beamx=xycenter(1) beamy=xycenter(2) type *, 'nx=',nx,' ny=',ny,' nbytes=',nbytes type *, 'beam center=',beamx, beamy IF (IPSYSTEM.EQ.'R-AXIS') THEN ENCODE(80,6739,TITLE) AQDATE,AQTIME PIX=0.105 ELSE ENCODE(80,6740,TITLE) AQDATE,AQTIME PIX=0.150 ENDIF 6739 FORMAT(' R-AXIS image of: ',2A20) 6740 FORMAT('MAR IPS image of: ',2A20) 72 FORMAT(' ',A80) WRITE(6,72)TITLE write(3,'(A)') '50 752 moveto ('//title//') show' ENCODE(80,6745,TITLE) ICCOUNTS(1),ICCOUNTS(2),EXPTIME 6745 FORMAT('IC COUNTS AT START, END:',2I9,' EXPOSURE TIME (sec)',I5) WRITE(6,72)TITLE write(3,'(A)') '50 732 moveto ('//title//') show' ENCODE(80,6741,TITLE) LAMBDA,DISTANCE,PHI_START,PHI_END 6741 FORMAT('LAMBDA: ',F11.4,' DISTANCE:',F9.2,' PHI:',2F9.2) WRITE (6,72)TITLE write(3,'(A)') '50 742 moveto ('//title//') show' c Save lambda, distance for calculating resolution: wl=lambda dist = distance IF (XM.GT.0) GOTO 87 TYPE *,'READING IMAGE. FIRST PASS: FINDING MIN, MAX DENSITY' DO 85 IY = 1,NY READ(2,155,end=860) (BLINE(II),II=1,nbytes) DO 80 IX=1,NX C TYPE *, ILINE(IX) X=ILINE(IX) call byteswapi2(x) IF (X.LT.0)X=X+R64K DMEAN=DMEAN+X IF (X.GT.DMAX) DMAX=X 80 IF (X.LT.DMIN) DMIN=X 85 CONTINUE nyt=ny GOTO 865 860 nyt=iy-1 865 CLOSE(UNIT=2) c use only number of record read in computing average 870 DMEAN=DMEAN/(NX*NYt) TYPE *,'DMIN=',DMIN,' DMAX=',DMAX,' DMEAN=',DMEAN ENCODE(80,6750,TITLE) DMIN,DMAX,DMEAN 6750 FORMAT('MARXER8 BIG V940107: Min/Max/mean density = ',2f8.0,f8.2) write(3,'(A)') '50 712 moveto ('//title//') show' OPEN (UNIT=2,FILE=INFILE,readonly,STATUS='OLD') C SKIP OVER HEADER FOR NOW: READ(2,155) (BLINE(II),II=1,124) X0=DMEAN/2 XM=DMEAN*3 87 SCALE=(XM-X0)/256 write(6,*)'xm,x0,scale=',xm,x0,scale if (j2.eq.0) then I1=400 I2=1599 J1=400 J2=1599 endif c CHECK THAT WINDOW LIMITS MAKE SENSE if (I1.gt.I2) then i=i1 i1=i2 i2=i endif if (I1.lt.1) then I1 = 1 write(6,*)'*****WARNING- WINDOW LIMIT RESET TO 1 *****' write(6,*)'pixels counted as 1 to ',ny,'*****' endif if (I1.ge.ny-1) I1=ny-2 if (I2.gt.ny) I2=ny if (J1.gt.J2) then i=j1 j1=j2 j2=i endif if (J1.lt.1) then J1 = 1 write(6,*)'*****WARNING- WINDOW LIMIT RESET TO 1 *****' write(6,*)'pixels counted as 1 to ',nx,'*****' endif if (J1.ge.nx-1) J1=nx-2 if (J2.gt.nx) J2=nx if ((I1.EQ.I2).OR.(I1.EQ.I2)) then WRITE(6,*)'IMAGE WINDOW HAS ZERO WIDTH OR HEIGHT!' WRITE(6,*) J1,J2,I1,I2 CLOSE(UNIT=2) GOTO 5 ENDIF in=1+i2-i1 jn=1+j2-j1 nmax=in if (jn.gt.nmax) nmax=jn ENCODE(80,6742,TITLE) X0,XM,I1,I2,J1,J2 6742 FORMAT('Density Limits:',2F8.1, . ' Selected Map Limits (X,Y)=',4I5) WRITE (6,72)TITLE write(3,'(A)') '50 762 moveto ('//title//') show' ENCODE(80,6743,TITLE) float(nmax)/20.,PIX*float(nmax)/20. c ENCODE(80,6743,TITLE) float(nmax)/19.5,PIX*float(nmax)/19.5 !for kodak printer 6743 FORMAT('Scale: 1 cm=',F8.2,' pixels =',f8.3, . ' mm on imaging plate') WRITE (6,72)TITLE write(3,'(A)') '50 722 moveto ('//title//') show' write(3,'(A)') '583 363.6 translate -90 rotate' write(3,'(A)') '0 4 moveto (TOP (view from source)) cshow' C WRITE POSTSCRIPT PRELUDE: WRITE(3,'(A)')'gsave initgraphics' WRITE(3,'(A)')' {1 exch sub} settransfer' nn=jn !nn is number of bytes in string- = # pix in record for 8-bit data J3=J2 c if (2*nn.lt.jn) then c nn=(jn+1)/2 c j3=j2+1 ! J3 MUST BE EVEN NUMBER FOR 4-BIT DATA. c OUTSTRING(J3)=HSTRING(1) ! PUT ASCII HEX CHAR IN LAST BYTE IN CASE ODD # PIXELS c endif WRITE(3,76) nn 76 format (' /imline ',i5,' string def') write(3,77) jn,in,nmax,-nmax,nmax 77 format (' /drawimage {',2i5,' 8 [',i5,' 0 0 ',i5,' 0 ',i5,']') WRITE(3,'(A)')' {currentfile imline readhexstring pop} image} def' c scale for unit square = 20 cm square c for 20 cm square, 20 * 72/2.54 = 566.9 write(3,'(A)')'18 75 translate 566.9 dup scale' c write(3,'(A)')'18 75 translate 552.8 dup scale' !for kodak printer. also change above. WRITE(3,'(A)')'drawimage' c note i1,i2,in refer to dimension across records. c MRC considers this Y, mars X c j1,j2,jn refer to dimension along one record. MRC X, mars Y. DO 88 IY = 1,I1-1 88 READ(2,155,end=6000) (BLINE(II),II=1,nbytes) DO 95 IY = I1,I2 READ(2,155,end=980) (BLINE(II),II=1,nbytes) DO 90 IX=J1,J3 call byteswapi2(ILINE(IX)) X=ILINE(IX) c IF (X.LT.0)X=X+R32K IF (X.LT.0)X=X+XM JX=INT((X-X0)/SCALE) IF (JX.GT.255) JX=255 IF (JX.LT.0) JX=0 hinib=int(jx/16) lonib= jx-16*hinib OUTSTRING(2*IX-1)=HSTRING(hinib+1) !1<=IX<=NX, OUTSTRING(1:2*NX) DIMENSIONED 90 OUTSTRING(2*IX)=HSTRING(lonib+1) 95 WRITE(3,156)(OUTSTRING(IX),IX=2*J1-1,2*J3) ! must be even number of bytes for 4-bit data, ALL ASCII HEX! 156 FORMAT(132A1) GOTO 990 980 WRITE(6,*)'UNEXPECTED END OF FILE READING IMAGE FILE! FILL W WHITE.' iyt=iy do 97 ix=j1,j3 97 outstring(ix)=hstring(1) do 98 iy=iyt,i2 98 WRITE(3,156)(OUTSTRING(IX),IX=J1,J3) 990 CLOSE(UNIT=2) C WRITE OUT REST OF POSTSCRIPT: WRITE(3,'(A)')'0 0 moveto 0 1 lineto 1 1 lineto 1 0 lineto 0 0 lineto' WRITE(3,'(A)')'0.000855 setlinewidth closepath clip' WRITE(3,'(A)')'{} settransfer stroke' z=1/float(nmax) c Scale for 1 unit = 1 pixel (now 1 unit = nmax pixels) WRITE(3,343) z,z 343 FORMAT(' 0 1 translate -90 rotate ',2F10.7 ' scale') c Should be -i1, -j1 translate, but denzo spots centered with 1 added. WRITE(3,344) 1-i1,1-j1 344 format(2i6,' translate ') WRITE(3,'(A)') '0 setlinewidth 0 setgray {} settransfer' IF (FNAME(:3).EQ.' ') GOTO 501 IF (BOXFLAG.EQ.' ') GOTO 501 IF (BOXFLAG.NE.'D') GOTO 501 c***************************** c***************************** WRITE(3,'(A)') '/cir { stroke r 0 360 arc stroke} def' WRITE(3,'(A)') '/sqr {moveto r dup neg dup rmoveto 2 mul dup 0 exch' WRITE(3,'(A)') ' rlineto dup 0 rlineto neg dup 0 exch rlineto' WRITE(3,'(A)') ' 0 rlineto closepath stroke} def' WRITE(3,'(A)') '/r 6 def' write(3,'(A)')'/Times-Roman findfont 6 scalefont setfont' write(3,*) '%1 1 .5 setrgbcolor' TYPE *,'READING DENZO OUTPUT FILE for spot locations' OPEN(UNIT=2,FILE=FNAME,STATUS='OLD') DO 345 I=1,5 345 READ (2,350)ASTRING C TYPE*, ASTRING 350 FORMAT (A40) N1=0 N2=0 C 355 READ (2,360)X 360 FORMAT (F7.4) !370 READ(2,101,END=500)H,K,L,J,X1,X2,X,S,X,X,Y,PL,X3 !101 FORMAT (I4,2I4,I2,F8.0,F8.1,F7.2,F6.1,F6.3,F7.1,F7.1,F6.3,F8.1) 370 READ(2,101,END=500)H,K,L,J,X1,X2,S,X,X,Y,PL,X3 101 FORMAT (I4,2I4,I2,F8.0,F8.1,7x,F6.1,F6.3,F7.1,F7.1,F6.3,F8.1) if (h.eq.999) goto 500 N1=N1+1 C IGNORE LOW-INTENSITY SPOTS C IF (X1.LT.100) GOTO 370 C Ignore offscale spots: if ((x.lt.i1).or.(x.gt.i2)) goto 370 if ((y.lt.j1).or.(y.gt.j2)) goto 370 IF (J.EQ.1) then C PUT CIRCLE AROUND PARTIAL WRITE(3,371) X,Y 371 FORMAT(2F7.1,' cir') else C PUT SQUARE AROUND FULLY RECORDED 380 WRITE(3,381) X,Y 381 FORMAT(2F7.1,' sqr') endif c label w hkl: write(3,383)x,y,h,k,l 383 format(2f7.1,' moveto (',3i3.0') show') GOTO 370 500 CLOSE (UNIT=2) write(3,'(A)')' stroke' C501 write(6,*)'Enter beam center x,y, wl, and distance; or 0,0,0,0', C & ' to not draw circles:' C write(6,*)'(beam center and distance in mm, wl in A.U.)' C read(5,*)beamx,beamy,wl,dist c501 write(6,*)'Enter beam center in mm: x,y; or 0,0 to not draw ', c & 'circles:' c read(5,*)beamx,beamy c if (beamx.eq.0) goto 505 c convert to pixels: c beamx=beamx/PIX c beamy=beamy/PIX c WRITE(6,*)' Assuming center is ',npix/2,npix/2,' pixels.' c BEAMX=npix/2 c BEAMY=npix/2 501 continue write(6,*)'Enter resolution in angstroms at which to draw circle:' write(6,*)'Enter 0 to draw no (more) circles' write(6,*)' wl=',wl write(3,*) '%0 1 1 setrgbcolor' 502 read(5,*)resol if (resol.eq.0) goto 505 c calc radius in mm x=wl/(2*resol) theta=ASIN(x) radius=dist*TAN(2*theta) c convert to pixels: radius=radius/PIX write(6,*)beamx,beamy,radius write(3,504) beamx,beamy,radius 504 format (' ',3f8.2,' 0 360 arc stroke') goto 502 505 WRITE(3,'(A)')' showpage grestore' CLOSE (UNIT=3) TYPE *,'NUMBER OF SPOTS REPORTED:',N1,' NUMBER USED:',N2 write(6,*)'LPR-ing marxer8.ps' C J=LIB$SPAWN('LPR MARXER8.PS') C J=LIB$SPAWN('DEL MARXER8.PS;*') c write(6,*)'purging marxer8.ps' c J=LIB$SPAWN('PURGE MARXER8.PS') c better not to purge, system may not have printed yet if que backed up c Instead, define lpr with /del so system wil delete c lpr :== print/del/param=:(data=post)/que=s07_ps_lps17b GOTO 5 6000 write(3,'(A)')' grestore 50 200 moveto (incomplete file!) show' close (unit=3) close (unit=2) WRITE(6,*) 'END OF FILE ENCOUNTERED BEFORE SELECTED AREA!' GOTO 505 END subroutine byteswapr4(value) real*4 value,exch1,exch2 character*1 exbyt1(4),exbyt2(4) equivalence (exch1,exbyt1(1)) equivalence (exch2,exbyt2(1)) exch1=value exbyt2(1)=exbyt1(2) exbyt2(2)=exbyt1(1) exbyt2(3)=exbyt1(4) exbyt2(4)=exbyt1(3) c value= exch2*4.0 !if run on vms-type machine (but file access needs fixing) value= exch2/4.0 !if run on unix machine return end c----------------------------------------- subroutine byteswapi4(value) integer*4 value,exch1,exch2 character*1 exbyt1(4),exbyt2(4) equivalence (exch1,exbyt1(1)) equivalence (exch2,exbyt2(1)) exch1=value exbyt2(1)=exbyt1(4) exbyt2(2)=exbyt1(3) exbyt2(3)=exbyt1(2) exbyt2(4)=exbyt1(1) value= exch2 return end c----------------------------------------- subroutine byteswapi2(value) integer*2 value,exch1,exch2 character*1 exbyt1(2),exbyt2(2) equivalence (exch1,exbyt1(1)) equivalence (exch2,exbyt2(1)) exch1=value exbyt2(1)=exbyt1(2) exbyt2(2)=exbyt1(1) value= exch2 return end