PROGRAM SPRAT C C A PROGRAM TO CALCULATE STOPPING POWERS AND RANGES BASED ON C THE OAK RIDGE SPAR CODE. C (ARMSTRONG & CHANDLER - ORNL 4869 [1973]) C C PROGRAM ADAPTED TO ACCEPT ENERGY STEP INSTEAD OF USING A C LOGARITHMIC SCALE , ALSO ALLOWS ALTERATIONS TO MEDIUM AND C ENERGY RANGE WHILE PROGRAM IS RUNNING . PROGRAM NOW GIVES C RANGE IN MG/CM2 AS WELL AS CM. C THE ORDER IN WHICH DATA IS INPUT WAS ALTERED BY REMOVING C THE REQUEST FOR A NAME FOR THE ION AND MOVING THE C TERMINATION REQUEST TO THE END OF THE MAIN ROUTINE. C THE ROUTINES 'EXEC' & 'SCHFL' ARE NOT AVAILABLE AT THIS C INSTALLATION AND HAVE BEEN COMMENTED OUT (HIGH LIGHTED BY A C ROW OF '*' ON EITHER SIDE) .THE ROUTINE RMPAR WAS REPLACED C BY A TERMINAL REQUEST FOR THE OUTPUT STREAM NUMBER. C (MANNING - 1981) C implicit double precision (a-h,o-z) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ DIMENSION NUCL(8),AMUN(4),XSN(4) DIMENSION MCHEM(10),ZMED(10),AMUA(10) C LOGICAL*1 MNAME(5) CHARACTER ANSWER*1,name*2 INTEGER OVDU C C MCHEM(I) - # ATOMS OF ELEMENT(I) IN MOLECULE OF MEDIUM C C ZMED(I) - ATOMIC NUMBER OF ELEMENT(I) C C AMUA(I) - MASS OF ONE ATOM OF ELEMENT(I) IN A.M.U. C DATA IVDU,OVDU/5,6/ LIST = 6 WRITE(OVDU,10) 10 FORMAT +('Do you want output redirected from screen file dedx_out + on stream 8 Y/N ?') READ(5,'(A1)') ANSWER IF (ANSWER.EQ.'Y'.OR.ANSWER.EQ.'y') THEN OPEN (8,FILE='dedx_out') LIST=8 ENDIF C IPAGE=1100B+LIST C IPAGE=12+LIST 220 WRITE(OVDU,1011) 1011 FORMAT('Figures are available for all chemical elements plus'/ 1'the compounds ISOBUTANE(IB) and MYLAR(MY). for any of the'/ 1'elements(standard isotopic composition) type its chemical'/ 1'symbol. for any other compound type (OT) ?') READ(IVDU,2000) name 2000 FORMAT(a2) * * CONVERT TO UPPER CASE * CALL UCASE(name,2) * PRESS = 760.0d0 NMED = 1 IF(NAME.EQ.'OT'.OR.NAME.EQ.'ot') GO TO 20 IF(NAME.NE.'MY'.AND.NAME.NE.'my') GO TO 77 NMED = 3 IANSW = 'N' AVDEN = 1.395d0 MCHEM(1) = 10 MCHEM(2) = 8 MCHEM(3) = 4 ZMED(1) = 6.0d0 ZMED(2) = 1.0d0 ZMED(3) = 8.0d0 AMUA(1) = 12.0d0 AMUA(2) = 1.0d0 AMUA(3) = 16.0d0 GO TO 30 77 IF(NAME.NE.'IB'.AND.NAME.NE.'ib') GO TO 78 NMED = 2 IANSW = 'Y' DNTP = 0.00259d0 WRITE(OVDU,*) 'Pressure (MMS)?' READ(IVDU,*) PRESS AVDEN = DNTP*PRESS/760.0 MCHEM(1) = 4 MCHEM(2) = 10 ZMED(1) = 6.0d0 ZMED(2) = 1.0d0 AMUA(1) = 12.0d0 AMUA(2) = 1.0d0 GO TO 30 78 CONTINUE CALL MELOOK(NAME,IANSW,MCHEM(1),AMUA(1),AVDEN,MZMED,FLAG,M) ZMED(1)=dble(MZMED) IF(AVDEN.NE.0.0d0) GO TO 101 WRITE(OVDU,*) 'Please type the density as this is unavailable?' READ(IVDU,*) AVDEN 101 IF(FLAG.EQ.0.0) GO TO 20 IF(IANSW.EQ.'N') GO TO 30 WRITE(OVDU,*) 'Pressure (MMS) ?' READ(IVDU,*) PRESS GO TO 30 C # OF DIFF ELEMENTS 20 WRITE(OVDU,1000) 1000 FORMAT('# of different elements in stopping medium ?') READ(IVDU,*)NMED C AVERAGE DENSITY WRITE(OVDU,1009) 1009 FORMAT('Medium gaseous [YES] or [NO] ?') READ(IVDU,2000)IANSW PRESS=760.0 IF(IANSW.EQ.'Y'.OR.IANSW.EQ.'y') GO TO 80 WRITE(OVDU,1001) 1001 FORMAT('Density of medium (GMS/CC) ?') READ(IVDU,*)AVDEN GO TO 90 80 WRITE(OVDU,1010) 1010 FORMAT('Density @ N.T.P.(GMS/CC), pressure (MMS) ?') READ(IVDU,*)DNTP,PRESS AVDEN=DNTP*PRESS/760.0d0 C C #ATOMS/MOL,Z,MASS FOR EACH ELEMENT C 90 DO 100 I=1,NMED WRITE(OVDU,1002)I 1002 FORMAT('ELEMENT #',I2/'# ATOMS IN MOLECULE, Z, MASS (AMU) ?') READ(IVDU,*)MCHEM(I),ZMED(I),AMUA(I) 100 CONTINUE C C AMUM - MASS OF MOLECULE IN A.M.U. C 30 AMUM=0.0 DO 110 I=1,NMED 110 AMUM=AMUM+MCHEM(I)*AMUA(I) C C DENM * 10**24 - MOLECULES/CC C DENM=AVDEN/(AMUM*1.660543) C C CALCULATE ELECTRONS/CC C AV. IONIZATION POTL C AV. CHARGE C SUMN=0.0d0 SUMNZ=0.0d0 SUMNZI=0.0d0 DO 120 I=1,NMED C EN - #ATOMS(I)/CC EN=MCHEM(I)*DENM C ENZ - #ELECTRONS(I)/CC ENZ=EN*ZMED(I) C ENZI - ENZ*ALOG(IONIZATION POTL(I)) ENZI=ENZ*ALGIP(ZMED(I)) SUMN=SUMN+EN SUMNZ=SUMNZ+ENZ 120 SUMNZI=SUMNZI+ENZI C ELNI - ALOG(MEAN IONIZATION POTL(MEV)) ELNI=SUMNZI/SUMNZ C AVIP - MEAN IONIZATION POTL(MEV) AVIP=EXP(ELNI) C EDEN * 10**24 - ELECTRONS/CC EDEN=SUMNZ C AVZ - MEAN CHARGE AVZ=SUMNZ/SUMN C C READ IN PARAMETERS OF BEAM TO BE ANALYZED C NNUC=1 C ********************************************************** C CALL SCHFL(NUCL,AMUN,XSN,NFND,NNUC) C IF(NFND(1).GE.0) GO TO 210 C WRITE(OVDU,1005) NUCL(1),NUCL(2) C1005 FORMAT(/2A2,' NOT FOUND, ENTER -VE MASS TO EXIT. MASS (AMU) ?') C READ(IVDU,*)AMUN(1) C ********************************************************** DO 625 I=1,5 625 MNAME(I) = ' ' 3010 WRITE(OVDU,1006) 1006 FORMAT('Please input the chemical symbol of the ion e.g 12C :-') READ(IVDU,4000) MNAME 4000 FORMAT(5A1) CALL GETM(M,MZMED,IFLAG,MNAME) ZPART=dble(MZMED) AMUN(1) = dble(M) IF(IFLAG.NE.1) GO TO 102 WRITE(OVDU,*)'Sorry unavailable please type MASS(AMU) & CHARGE?' READ(IVDU,*) AMUN(1),ZPART C IF(AMUN(1).LE.0.0) GO TO 999 102 XSN(1)=0.0 210 AMASS=AMUN(1)+XSN(1)/931481.2d0 EMASS=AMASS*931.4812d0 C C PRINT OUT DETAILS OF MEDIUM C C ********************************************************* C CALL EXEC(3,IPAGE,-1) C ********************************************************* IF(IANSW.NE.'Y') GO TO 225 WRITE(LIST,3003)NAME,AVDEN,PRESS 3003 FORMAT(15X,'MEDIUM IS ',A2//G15.6,' GMS/CC , ',G15.6,' MM'/) GO TO 3007 225 WRITE(LIST,3006) NAME,AVDEN 3006 FORMAT(15X,'MEDIUM IS ',A2//10X,G15.6,' GMS/CC') 3007 CONTINUE WRITE(LIST,3005) (MCHEM(I),ZMED(I),AMUA(I),I=1,NMED) 3005 FORMAT(' MOLECULE:'/(1X,I4,' ATOMS OF Z=',F10.6,' ,MASS=',F10.5, 1 ' AMU')) WRITE(LIST,3004)AVIP,EDEN,AVZ 3004 FORMAT(/'AVE IONIZATION POTL IS ',G15.6/'ELECTRON DENSITY*10**24' 1 ,G15.6/'AVERAGE CHARGE ',G15.6/) WRITE(LIST,3000)AMASS,ZPART 3000 FORMAT(' ION:'/' ',F7.3,' AMU , Z =',F7.3/) 3008 WRITE(OVDU,1007) 1007 FORMAT('emin,emax(MEV),estep ?') READ(IVDU,*) EMIN,EMAX,ESTEP IF(EMIN.EQ.0.) EMIN = 0.01 IF(EMAX.LT.EMIN)EMAX=EMIN IF(ESTEP.LE.0) ESTEP = 1 WRITE(OVDU,*) 'Do you require the range values?' READ(IVDU,2000) JANSW IF(JANSW.EQ.'Y'.OR.JANSW.EQ.'y') GOTO 18 IF(IANSW.EQ.'N'.OR.IANSW.EQ.'n') GOTO 19 WRITE(LIST,*) 'ENERGY(MEV) DE/DX(MEV/CM) DE/DX(MEV.CM2/MG)' GOTO 60 19 WRITE(LIST,*) 'ENERGY(MEV) DE/DX(MEV/MM) DE/DX(MEV.CM2/MG)' GOTO 60 18 IF(IANSW.EQ.'N'.OR.IANSW.EQ.'n') GO TO 50 WRITE(LIST,3002) 3002 FORMAT(' ENERGY(MEV) DE/DX(MEV/CM) DE/DX(MEV.CM2/MG)', 3 ' RANGE(MM) RANGE(MG/CM2)'/) GO TO 60 50 WRITE(LIST,3009) 3009 FORMAT(' ENERGY(MEV) DE/DX(MEV/MM) DE/DX(MEV.CM2/MG)', 3 ' RANGE(MM) RANGE(MG/CM2)'/) 60 EPART=EMIN ETEMP = 0.0d0 DX = 0.5d0/DEDX(EMASS,EPART,ZPART) 226 SP=DEDX(EMASS,EPART,ZPART) SPK=(SP*0.001)/AVDEN IF(JANSW.EQ.'Y'.OR.JANSW.EQ.'y') GOTO 16 IF(IANSW.EQ.'Y'.OR.IANSW.EQ.'y') GOTO 17 WRITE(LIST,3011) EPART,SP/10,SPK GOTO 25 17 WRITE(LIST,3011) EPART,SP,SPK GOTO 25 16 R=RANGE(DX,EMASS,EPART,ZPART,SP) RCM = R*AVDEN*1000.0d0 IF(IANSW.EQ.'Y'.OR.IANSW.EQ.'y') GO TO 15 WRITE(LIST,3001) EPART,SP/10,SPK,R*10.d0,RCM GOTO 25 15 WRITE(LIST,3001)EPART,SP,SPK,R,RCM 3001 FORMAT(5G15.6) 3011 FORMAT(3G15.6) 25 DX=R/50 ETEMP = ETEMP+1.0d0 EPART = EMIN +(ETEMP*ESTEP) IF(EPART.LE.EMAX) GO TO 226 230 CONTINUE 200 WRITE(OVDU,1003) 1003 FORMAT(/'/ TO EXIT,M= NEW MEDIUM,L= NEW LIMITS,I= NEW ION,', 1 'R= RANGE VALUES ?') READ(IVDU,1004) NUCL(1),NUCL(2) 1004 FORMAT(2A2) IF(NUCL(1).EQ.'M'.OR.NUCL(1).EQ.'m') GO TO 220 IF(NUCL(1).EQ.'L'.OR.NUCL(1).EQ.'l') GO TO 3008 IF(NUCL(1).EQ.'I'.OR.NUCL(1).EQ.'i') GO TO 3010 IF(NUCL(1).NE.'R'.AND.NUCL(1).NE.'r') GO TO 999 JANSW = 'Y' GOTO 18 999 STOP END FUNCTION DEDX(EMASS,EPART,ZPART) C C ONLY VALID FOR BETA > 0.0046*Z**(1/3) C implicit double precision (a-h,o-z) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ * DOUBLE PRECISION BETA2,DB2 DATA ZP/1.0d0/ DB2=BETA2(EPART,EMASS) BETA=sqrt(DB2) PE=BTOEP(DB2) GOG=ZEFF(ZPART,BETA)/ZEFF(ZP,BETA) GOG=GOG*GOG DEDX=GOG*DEDXP(PE,DB2) RETURN END FUNCTION BTOEP(DBSQ) C C TO CALCULATE THE ENERGY OF A PROTON GIVEN BETA**2, WHERE C BETA IS THE SPEED OF THE PROTON RELATIVE TO THE SPEED OF LIGHT C C INPUT:- DBSQ - BETA**2 (DOUBLE PRECISION) C implicit double precision (a-h,o-z) * DOUBLE PRECISION DBSQ,D DATA EMP/938.2592d0/ D=sqrt(1D0/(1D0-DBSQ))-1D0 BTOEP=EMP*D RETURN END FUNCTION RANGE(DX,EMASS,EPART,ZPART,SP) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ implicit double precision (a-h,o-z) * DOUBLE PRECISION ED SPE=SP ED=EPART R=0.0d0 ES=EPART 100 ED=ED-DE(DX,EMASS,ES,ZPART,SPE) E=ED IF(E.LE.0.0000001d0) GO TO 200 R=R+DX ES=E SPE=DEDX(EMASS,ES,ZPART) GO TO 100 200 RANGE=R+ES/SPE RETURN END FUNCTION DE(DX,EMASS,EPART,ZPART,SP) C C TO CALCULATE ENERGY LOST OVER DX CMS ASSUMING A QUADRATIC C RELATIONSHIP BETWEEN E & X. C implicit double precision (a-h,o-z) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ DELTAE=SP*DX ENEW=EPART-DELTAE IF(ENEW.LE.0.0d0) GO TO 200 SPGES=DEDX(EMASS,ENEW,ZPART) DE=DX*(0.75d0*SP+(0.25d0*SPGES/SP)*SPGES) RETURN 200 DE=EPART RETURN END FUNCTION DEDXP(ENRGY,DB2) C C ONLY VALID FOR BETA > 0.0046 : 10 KEV PROTONS C implicit double precision (a-h,o-z) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ * DOUBLE PRECISION D,DB2,DSP DATA EMASS/938.2592d0/,Z/1.0d0/ BETA=sqrt(DB2) C C SP=(EDEN*12.56637*532.281)/(0.5110041*2.567018) C C *10**24 4*PI *10**-40 ME*C**2 *10**-12 C CM**-3 E**4 ESU MEV (ERG/MEV)**2 C C DSP=SP*1.0D-4 C DSP=(EDEN*5.099147D-1) ZE=ZEFF(Z,BETA) DSP=DSP*(ZE*ZE)/DB2 D=log(1.022008D0*DB2/(1.0D0-DB2))-DB2 DELTA=DEFF(ENRGY) COZ=SHELL(ENRGY) DSP=DSP*(D-ELNI-COZ-DELTA/2.0D0) DEDXP=DSP RETURN END FUNCTION SHELL(E) C C TO CALCULATE THE SHELL CORRECTION TERM IN STOPPING POWER C CALCULATIONS(SPAR-ARMSTRONG & CHANDLER-ORNL-4869 [1973]) C C INPUT - E ENERGY IN MEV C AVIP MEAN IONIZATION POTENTIAL OF STOPPING MEDIUM (MEV) C AVZ MEAN ATOMIC NUMBER OF STOPPING MEDIUM C implicit double precision (a-h,o-z) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ * DOUBLE PRECISION BETA2 * DOUBLE PRECISION GNU2 * DOUBLE PRECISION F1,F2 DATA A1/0.422377d-6/,A2/3.858019d-9/, 1 B1/0.304043d-7/,B2/-0.1667989d-9/, 2 C1/-0.38106d-9/,C2/0.157955d-11/ DATA EMP,ZAL,ZH2O/938.2592d0,12.95d0,3.34d0/ DATA P1,P2,P3,P4,P5/4.774248d-4,1.143478d-3,-5.63392d-2, 1 4.763953d-1,4.844536d-1/ DATA W1,W2,W3,W4,W5,W6,W7,W8/-1.819954d-6,-2.232760d-5, 1 1.219912d-4,1.837873d-3,-4.457574d-3,-6.837103d-2, 2 5.266586d-1,3.743715d-1/ IF(E.LT.8.0) GO TO 10 GNU2=(E/EMP)*(E/EMP+2.0d0) F1= A1/GNU2 + B1/(GNU2*GNU2) + C1/(GNU2**3) F2= A2/GNU2 + B2/(GNU2*GNU2) + C2/(GNU2**3) SHELL= AVIP*AVIP*(F1 + F2*AVIP)/AVZ GO TO 999 10 B2=BETA2(E,EMP) SHELL=LOG(1.0220082d0) + LOG(B2) - LOG(AVIP) X=18769.0d0*B2/AVZ XLOG=LOG(X) XLOG2=XLOG*XLOG XLOG3=XLOG2*XLOG XLOG4=XLOG2*XLOG2 IF(AVZ.GT.ZAL) GO TO 20 C Z < ZAL XLOG5=XLOG2*XLOG3 XLOG6=XLOG3*XLOG3 XLOG7=XLOG3*XLOG4 XL1=W1*XLOG7 + W2*XLOG6 + W3*XLOG5 + W4*XLOG4 + W5*XLOG3 1 + W6*XLOG2 + W7*XLOG + W8 XL1=EXP(XL1) IF(AVZ.GT.ZH2O) GO TO 20 C Z < ZH2O XL=XL1 GO TO 30 C Z > ZH2O 20 XL=P1*XLOG4 + P2*XLOG3 + P3*XLOG2 + P4*XLOG + P5 XL=EXP(XL) IF(AVZ.GT.ZAL) GO TO 30 C ZH2O < Z < ZAL XL2=XL XL=XL1 + (AVZ-ZH2O)/(ZAL-ZH2O) * (XL2-XL1) C C PUT IN L(X) C 30 SHELL=SHELL - XL 999 RETURN END FUNCTION DEFF(E) C C TO CALCULATE THE DENSITY-EFFECT CORRECTION TERM IN STOPPING C POWER CALCULATIONS (SPAR-ARMSTRONG & CHANDLER-ORNL-4869 [1973]) C C INPUT - E ENERGY IN MEV C ELNI ALOG(MEAN IONIZATION POTENTIAL OF MEDIUM (MEV) ) C EDEN ELECTRON DENSITY OF STOPPING MEDIUM C implicit double precision (a-h,o-z) COMMON ISECT(128),EDEN,ELNI,AVIP,AVZ * DOUBLE PRECISION DEL DATA EMP/938.2592d0/ DEL=E/EMP DEL=log(DEL) + log(DEL+2.0d0) DEL=log(1.378D-9*EDEN) + DEL - 2.0D0*ELNI - 1.0D0 DEFF=DEL IF(DEFF.LT.0.0d0) DEFF=0.0d0 RETURN END DOUBLE PRECISION FUNCTION BETA2(E,EM) C C TO CALCULATE BETA**2 WHERE BETA IS THE SPEED OF A PARTICLE C RELATIVE TO THE SPEED OF LIGHT C C INPUT - E ENERGY IN MEV C EM REST MASS IN MEV C implicit double precision (a-h,o-z) * DOUBLE PRECISION R R=E/EM + 1D0 BETA2=1D0-(1D0/(R*R)) RETURN END FUNCTION ZEFF(Z,BETA) C C TO CALCULATE THE EFFECTIVE CHARGE OF A PARTICLE IN C STOPPING POWER CALCULATIONS(SPAR-ARMSTRONG & CHANDLER-ORNL- C 4869 [1973]) C C INPUT - Z NOMINAL CHARGE OF NUCLEUS C BETA SPEED OF THE PARTICLE RELATIVE TO THE SPEED OF LIGHT C implicit double precision (a-h,o-z) ZP=Z**0.666667d0 IF(BETA.LE.(0.07d0*ZP))GO TO 100 ZEFF=Z RETURN 100 ZEFF= Z * (1.0d0 - EXP(-125.0d0*BETA/ZP)) RETURN END FUNCTION ALGIP(Z) C C TO CALCULATE ALOG(IONIZATION POTENTIAL) FOR AN ELEMENT C OF ATOMIC NUMBER Z C C N.B. IONIZATION POTL IN MEV !!!!!! C implicit double precision (a-h,o-z) DIMENSION POT(13) DATA POT/18.7d0,42.0d0,39.0d0,60.0d0,68.0d0,78.0d0,99.5d0, & 98.5d0,117.0d0,140.0d0, & 150.0d0,157.0d0,163.0d0/ IZ=(Z+0.05d0) IF(IZ.GT.13) GO TO 100 POTL=POT(IZ) GO TO 200 100 POTL=9.76d0*Z + 58.8d0/(Z**0.19d0) 200 ALGIP=LOG(POTL*1.0d-6) RETURN END SUBROUTINE MELOOK(NAME,IANSW,NO,AMASS,AVDEN,MZMED,FLAG,M) double precision amass,avden,flag character*2 name,nname INTEGER IOSTR IOSTR = 11 FLAG = 1.0d0 c CALL CONECT('/home/oef/Research/anprogs/lib/spdata.dat',FILE) OPEN(IOSTR,STATUS='OLD', +FILE='/home/oef/Research/anprogs/lib/dedx.dat') 10 READ(IOSTR,30,END=100) IANSW,NO,NNAME,MZMED,AMAS_r,AVDE_r,M 30 FORMAT(A2,I1,A2,I7,F14.6,F14.8,I7) IF(NAME.NE.NNAME) GO TO 10 GO TO 999 100 FLAG = 0.0d0 999 close(IOSTR) amass=dble(amas_r) avden=dble(avde_r) RETURN END C *SUBROUITNE GETM TO INTERPRET A CHEMICAL SYMBOL IN 12C SYMBOLIC* C * FORMAT AND RETURN ATOMIC MASS AND ATOMIC NUMBER * C * INPUT.........MNAME(5),READ IN AS A5 FORMAT * C * OUTPUT........IA ATOMIC MASS * C * IZ ATOMIC NUMBER * C * IFLAG (0) O.K. * C * (1) MASS NOT FOUND IN NUCNAM* C C C C SUBROUTINE GETM(IA,IZ,IFLAG,MNAME) LOGICAL*1 NUM(3),LET(2),NOS(10),MNAME(5),BLANK INTEGER*2 ATOM,NUCNAM(110) CHARACTER FORMAT(3)*4,BUFF*24 EQUIVALENCE (ATOM,LET(1)) DATA NUCNAM/' N','H ','HE', +'LI','BE','B ','C ','N ','O ','F ','NE', +'NA','MG','AL','SI','P ','S ','CL','AR', +'K ','CA','SC','TI','V ','CR','MN','FE','CO', +'NI','CU','ZN','GA','GE','AS','SE','BR','KR', +'RB','SR','Y ','ZR','NB','MO','TC','RU','RH', +'PD','AG','CD','IN','SN','SB','TE','I ','XE', +'CS','BA','LA','CE','PR','ND','PM','SM', +'EU','GD','TB','DY','HO','ER','TM','YB', +'LU','HF','TA','W ','RE','OS','IR','PT', +'AU','HG','TL','PB','BI','PO','AT','RN', +'FR','RA','AC','TH','PA','U ','NP','PU', +'AM','CM','BK','CF','ES','FM','MD','NO', +'LR','RF','HA','NH','NS','UO','UE'/ DATA NOS/'0','1','2','3','4','5','6','7','8','9'/ DATA BLANK/' '/ DATA FORMAT /'(I1)','(I2)','(I3)'/ C C *RECOGNISING A AND Z IN MNAME* C CALL UCASE(MNAME,5) 10 DO 40 I=1,5 DO 30 K=1,10 IF (MNAME(I).EQ.NOS(K)) GOTO 40 30 CONTINUE GOTO 50 40 CONTINUE 50 DO 60 I2=1,2 LET(I2)=BLANK 60 CONTINUE M=I+1 J1=0 DO 70 L=I,M J1=J1+1 LET(J1)=MNAME(L) 70 CONTINUE N=I-1 J1=0 DO 80 L=1,N J1=J1+1 NUM(J1)=MNAME(L) 80 CONTINUE LENFLD=J1 C CALL DEFBUF(13,LENFLD,NUM) C READ(13,*)IA WRITE(BUFF,'(3A1)') NUM READ(BUFF,FORMAT(LENFLD)) IA C C *LOOKING FOR NUCLEI SYMBOL IN PERIODIC TABLE* C DO 90 I1=1,107 IF (ATOM.EQ.NUCNAM(I1)) GOTO 110 90 CONTINUE IFLAG=1 RETURN 110 IZ=I1-1 RETURN END subroutine ucase(string,in) * * subroutine to lowercase a string for a unix filename * character string*(*) * in=index(string,' ')-1 do 10 i=1,in ic=ichar(string(i:i)) if (ic.ge.97 .and. ic.le.122) string(i:i)=char(ic-32) 10 continue return end