From 5ba963c1ed94b55bed0efda768835237babf0c77 Mon Sep 17 00:00:00 2001 From: nemu Date: Fri, 22 Nov 2013 15:47:50 +0000 Subject: [PATCH] Added datmak7l to obsolete directory; added ranlux source to trimsp7l_tp.F. --- trimsp/obsolete/Makefile | 2 +- trimsp/obsolete/datmak7l_tp.F | 1007 +++++++++++++++++++++++++++++++++ trimsp/obsolete/trimsp7l_tp.F | 295 ++++++++++ 3 files changed, 1303 insertions(+), 1 deletion(-) create mode 100644 trimsp/obsolete/datmak7l_tp.F diff --git a/trimsp/obsolete/Makefile b/trimsp/obsolete/Makefile index 106f77e..731cc49 100644 --- a/trimsp/obsolete/Makefile +++ b/trimsp/obsolete/Makefile @@ -21,7 +21,7 @@ trimsp7l_tp.o : trimsp7l_tp.F $(FC) $(OPS) $< trimsp7l_tp : trimsp7l_tp.o - $(FC) -o $(HOME)/bin/$@ $< $(CERNLIB_DIR)/libmathlib.a + $(FC) -o $(HOME)/bin/$@ $< datmak7l_tp.o : datmak7l_tp.F $(FC) $(OPS) $< diff --git a/trimsp/obsolete/datmak7l_tp.F b/trimsp/obsolete/datmak7l_tp.F new file mode 100644 index 0000000..d1ff15c --- /dev/null +++ b/trimsp/obsolete/datmak7l_tp.F @@ -0,0 +1,1007 @@ +c + PROGRAM DATMAKER +c +c------------------------------------------------------------------------------ +c dieses Programm erstellt die Eingabedateien, die mit dem Programm +c TRIMSP4L gelesen werden. +c dies ist die Version DATMAK7L(-test) +c +c Festlegen der für TRIMSP benötigten Modellparameter +c Bezeichnung nach Eckstein's Datei TRVMC95.text +c ef,esb,sheath,erc,rd,ca,kk0,kk0r,kdee1,kdee2,ipot,ipotr,irl +c +c Projektilparameter: +c =================== +c +c zproj = Ordnungszahl des Projektils == fuer myonen = 1 +c mproj = Massenzahl des Projektils == fuer Myonen = 0.11 +c alphaproj = Winkel zwischen Projektil und Targetnormale == 0 ist +c senkrecht +c eproj = Projektilenergie in eV +c esig = Breite der gaussfoermigen Energieverteilung +c thick(n) = Dicke eines Layers in Angstroem +c rho(n) = Dichte eines Layers in g/cm^3 +c +c Simulationsparameter: +c ===================== +c +c ri = Zufallszahleninitialisierung (muss ungerade sein) +c ri2 = Zufallszahleninitialisierung, mit der die Energie gewuerfelt wird +c x0 = von wo aus werden Projektile implantiert == 0 ist Probenoberflaeche +c cw = Dickeninterval in Angstroem +c Festlegen der Target-Elemente m(max. 5) in den jeweiligen Layern l(max. 7) +c p1 = Ordnungszahl-->ordz(n) +c p2 = Massenzahl (amu)-->mass(n) +c p3 = Bindungsenergie im Festkoerper (eV)-->elas(n) +c p4 = Dichte des Elementes (g/cm^3)-->rho(n) +c p(1-5) = Stoppingpowerkoeffizienten fuer H in Materie nach +c Andersen-Ziegler-->spower(k,m,l) +c conz(m,l) --> Konzentration eines Elementes in einem Layer +c conzlayer(n) --> Gesamasskonzentration im Layer == muss 1.0 sein +c die Parameter SBE(5) sind auf 30.0 eV und BE(5) auf 0 gesetzt (siehe TRVMC95.txt) +c die Parameter CK(3) sind auf 1 gesetzt (siehe TRVMC95.txt) +c +c April 1999: das Programm berechnet jetzt aus den Konzentrationen die +c stöchiometrische Dichte eines Layers. Wichtig: die Dichten fuer +c H,He,Ne,O,Ar,Kr,Xe sind die Dichten am Tripelpunkt fuer die +c fluessige Phase, die Dichten fuer F,Cl,Rn,Fr,Ra sind nicht bekannt +c und deshalb auf 0.1 gesetzt. +c +c Mai 1999 : das Programm wurde erweitert. Jetzt koennen in einer Schleife +c mehrere Eingabefiles fuer TRIMSP (ab Version TRIMSPPc4i) mit +c verschiedener Startenergie und verschiedenen Kanalgroessen fuer +c die Range Profile erstellt werden. Ausserem weden die erstellten +c Files auf Abfrage in eine fuer ab TRIMSPpc4h geeignete +c Batch-Datei weggeschrieben. +c Juni 1999 : die Stopkoeffizienten aus der ICRU49 Tabelle werden verwendet. +c Datei :stopicru.dat +c Dez. 1999 : Dichten fuer N2,O2, Ar,Xe,Kr,Ne nun fuer die feste Phase +c : Batchtype x eingefuegt +c ab Version datmak4k +c +c ! Eingabedateien, die mit datmak4k erzeugt werden, koennen von aelteren Versionen (vor TrimSpp4k) +c ! nicht gelesen werden. +c +c Dez. 1999 : neu Variable +c alphaprojsig (Breite einer Gaussverteilung der Einfallwinkel +c ri3 = seed fuer Zufallszahlengenerator fuer Einfallwinkel +c Batchtype a,b eingefuegt aber noch nicht implementiert +c Jun. 2000: nun sieben layer +c +c major revision: +c +c 19-Jun-2002: TP, implicit none, program comments at the beginning, +c f77 standard, english output, changed formatting +c change extension to .F in order to invole the pre- +c compiler to write either a DOS bat file or a bash +c shell script. +c no longer prompt for random generator seed because +c ranlux is used in trimsp7l_tp. +c +c 14-Jan-2003: TP, replaced format '(I6)' of reading the number of projectiles nproj +c by *. +c +c-------------------------------------------------------------------------------------------- +c +c check OS +c +#if defined( _WIN32 ) +#define OS_WIN +#else +#define OS_UNIX +#endif +c + IMPLICIT none + + INTEGER i,k,l,m,lnum + INTEGER nproj,n,lmax + INTEGER kk0,kk0r,kdee1,kdee2,ipot,ipotr,irl + + REAL*4 zproj,mproj,eproj,esig,alphaproj,alphaprojsig + REAL*4 ri,ri2,ri3,cw,x0 + REAL*4 ef,esb,sheath,erc,rd,ca + REAL*4 p2,p3,p4,p(5) + REAL*4 ordz(5,7),mass(5,7),elas(5,7),spower(5,5,7),conz(5,7) + REAL*4 conzlayer(7),thick(7),rho(7) + REAL*4 ck(7),sbe(5,7),be(5,7) + + CHARACTER errcom*17,en*2,f*1,ausgabe*12,inpnam*4,parchar*3, + 1 inpext*4 + CHARACTER batchtype*1,name*8 + DATA zproj,mproj /1.0,0.11/ + DATA alphaproj,alphaprojsig /0.0,0.0/ + DATA eproj,esig /1000.0,0.0/ + DATA ef,esb,sheath,erc /0.50,0.00,0.00,0.00/ + DATA rd,ca /50.00,1.00/ + DATA kk0,kk0r,kdee1,kdee2,ipot,ipotr,irl /2,2,4,3,2,1,0/ + DATA ri,x0,cw /78741.00,0.00,30.0/ + DATA nproj /1000/ + DATA lmax /7/ + + DATA errcom /'Wrong Input !'/ + DATA inpext /'.inp'/ + DATA batchtype /'e'/ +c + WRITE(*,*)' ----------------------------------------' + WRITE(*,*)' | Program D A T M A K, Version 7L |' + WRITE(*,*)' ----------------------------------------' + WRITE(*,*) + WRITE(*,*)' Program to create the input file for' + WRITE(*,*)' TRIM.SP (W. Eckstein, IPP Garching).' + write(*,*) + WRITE(*,*)' The general name for the output (must be '// + 1 'character*4, A4) will be used' + WRITE(*,*)' in all output file names. In case of multiple'// + 1 ' energies to be simulated (Batchtype e),' + WRITE(*,*)' the energy will be appended as character*3 (A3),'// + 1 ' for example, 30keV=300, 0.5keV=005.' + WRITE(*,*)' Otherwise, the character*3 file name specifier'// + 1 ' must be entered manually.' + WRITE(*,*)' Between general file name and character*3'// + 1 ' specifier the batch type will be inserted' + WRITE(*,*)' as character*1.' + WRITE(*,*)' In batch mode only one of the paramters energy,'// + 1 'energy width, number of muons or layer thickness' + WRITE(*,*)' can be changed.' + WRITE(*,*) + WRITE(*,*)' Enter the general file name (A4)' + READ(5,'(A4)')inpnam +c +800 WRITE(*,*)' Batchtype ?' + WRITE(*,*)' E = energy (default)' + WRITE(*,*)' S = sigma energy' + WRITE(*,*)' N = number of particles' + WRITE(*,*)' D = layer thickness' + WRITE(*,*)' X = energy and sigma energy' +c WRITE(*,*)' A = alpha and alphasigma' +c WRITE(*,*)' B = energy, sigma energy, alpha, alphasigma' +c + WRITE(*,*)'Change batch type (y/n)? ' + READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + batchtype='e' + GOTO 806 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN +805 WRITE(*,*)' Enter batch type (A1) ' + READ(*,'(A1)')batchtype + IF(batchtype.EQ.'e'.OR.batchtype.EQ.'E')THEN + batchtype='e' + ELSEIF(batchtype.EQ.'s'.OR.batchtype.EQ.'S')THEN + batchtype='s' + ELSEIF(batchtype.EQ.'n'.OR.batchtype.EQ.'N')THEN + batchtype='n' + ELSEIF(batchtype.EQ.'d'.OR.batchtype.EQ.'D')THEN + batchtype='d' + ELSEIF(batchtype.EQ.'x'.OR.batchtype.EQ.'X')THEN + batchtype='x' +c ELSEIF(batchtype.EQ.'a'.OR.batchtype.EQ.'A')THEN +c batchtype='a' +c ELSEIF(batchtype.EQ.'b'.OR.batchtype.EQ.'B')THEN +c batchtype='b' + ELSE + WRITE(*,*)errcom + GOTO 805 + ENDIF + ELSE + WRITE(*,*)errcom + GOTO 800 + ENDIF +c +806 WRITE(*,*)' Projectile parameter ' + WRITE(*,'(1x,A3,F7.2,1x,A4,F7.2)')'Z: ',zproj,'M: ',mproj + WRITE(*,'(1x,A3,F7.2,1x,A6,F7.2)')'E: ',eproj, + + 'Esig: ',esig + WRITE(*,'(1x,A8,F7.2,1x,A10,F7.2)')'alpha: ',alphaproj, + + 'alphasig: ',alphaprojsig + WRITE(*,*)'Change projectile ? ' +1001 READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 1000 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' Charge number (Myon = 1)' +810 READ(5,501,ERR=810)zproj + IF(zproj.EQ.0.0) zproj=1. + WRITE(*,*)' Projectile mass (Myon=0.113)' +811 READ(5,501,ERR=811)mproj + IF(mproj.EQ.0.0) mproj=0.113 + ELSE + WRITE(*,*)errcom + GOTO 1001 + ENDIF +1000 WRITE(*,*)'Change projectile energy ? ' + READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 10011 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' Projectile energy in eV ' +814 READ(5,501,ERR=814)eproj + IF(eproj.EQ.0.0) eproj=1000. + ELSE + WRITE(*,*)errcom + GOTO 1000 + ENDIF +10011 WRITE(*,*)'Change sigma ? ' + READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 10012 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' Sigma in eV ' +815 READ(5,501,ERR=815)esig + esig=ABS(esig) + ELSE + WRITE(*,*)errcom + GOTO 10011 + ENDIF +10012 WRITE(*,*)'Change angle of projectile ? ' + READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 10013 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' Projectile angle (perpendicular=0)' +813 READ(5,501,ERR=813)alphaproj + alphaproj=ABS(alphaproj) + ELSE + WRITE(*,*)errcom + GOTO 10012 + ENDIF +10013 WRITE(*,*)'Change sigma ? ' + READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 10014 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' Sigma projectile angle ' +8130 READ(5,501,ERR=8130)alphaprojsig + alphaprojsig=ABS(alphaprojsig) + ELSE + WRITE(*,*)errcom + GOTO 10013 + ENDIF +c +10014 WRITE(*,*)' Other projectile parameters + + (see TRVMC95.txt for details)' + WRITE(*,'(1x,A3,F7.2,1x,A4,F7.2,1x,A7,F7.2,1x,A4,F7.2)') + + 'EF ',ef,'ESB ',esb,'SHEATH ',sheath,'ERC ',erc + WRITE(*,'(1x,A3,F7.2,1x,A3,F7.2)')'RD ',rd,'CA ',ca + WRITE(*,*)'Change ?' +1003 READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 1002 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' EF :' +820 READ(5,501,ERR=820)ef + WRITE(*,*)' ESB :' +821 READ(5,501,ERR=821)esb + WRITE(*,*)' SHEATH :' +822 READ(5,501,ERR=822)sheath + WRITE(*,*)' ERC :' +823 READ(5,501,ERR=823)erc + WRITE(*,*)' RD :' +824 READ(5,501,ERR=824)rd + WRITE(*,*)' CA :' +825 READ(5,501,ERR=825)ca +501 FORMAT(F7.0) + ELSE + WRITE(*,*)errcom + GOTO 1003 + ENDIF +1002 WRITE(*,*)' Elastic collision parameters + + (see TRVMC95.txt for details)' + WRITE(*,'(1x,A4,I1,1x,A5,I1,1x,A6,I1,1x,A6,I1)') + + 'KK0 ',KK0,'KK0R ',KK0R,'KDEE1 ',KDEE1,'KDEE2 ',KDEE2 + WRITE(*,'(1x,A5,I1,1x,A6,I1,1x,A4,I1)') + + 'IPOT ',IPOT,'IPOTR ',IPOTR,'IRL ',IRL + WRITE(*,*)'Change ?' +1005 READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 1004 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' KK0 ' +830 READ(5,502,ERR=830)KK0 + IF(KK0.LT.0.OR.KK0.GT.4)THEN + WRITE(*,*)' must be in [0,4].' + GOTO 830 + ENDIF + WRITE(*,*)' KK0R ' +831 READ(5,502,ERR=831)KK0R + IF(KK0R.LT.0.OR.KK0R.GT.4)THEN + WRITE(*,*)' must be in [0,4]. ' + GOTO 831 + ENDIF + WRITE(*,*)' KDEE1 ' +832 READ(5,502,ERR=832)KDEE1 + IF(KDEE1.LT.1.OR.KDEE1.GT.5)THEN + WRITE(*,*)' must be in [1,5].' + GOTO 832 + ENDIF + WRITE(*,*)' KDEE2 ' +833 READ(5,502,ERR=833)KDEE2 + IF(KDEE2.LT.1.OR.KDEE2.GT.3)THEN + WRITE(*,*)' must be in [1,3].' + GOTO 833 + ENDIF + WRITE(*,*)' IPOT ' +834 READ(5,502,ERR=834)IPOT + IF(IPOT.LT.1.OR.IPOT.GT.3)THEN + WRITE(*,*)' must be in [1,3].' + GOTO 834 + ENDIF + WRITE(*,*)' IPOTR ' +835 READ(5,502,ERR=835)IPOTR + IF(IPOTR.LT.1.OR.IPOTR.GT.3)THEN + WRITE(*,*)' must be in [1,3].' + GOTO 835 + ENDIF + WRITE(*,*)' IRL ' +836 READ(5,502,ERR=836) IRL + IF(IRL.LT.0.OR.IRL.GT.1)THEN + WRITE(*,*)' must be in [0,1]' + GOTO 836 + ENDIF +502 FORMAT(I1) + ELSE + WRITE(*,*)errcom + GOTO 1005 + ENDIF +c +1004 WRITE(*,*)' Parameter for the simulation:' + WRITE(*,'(1x,A18,I6)')'Number of projectiles ',nproj + WRITE(*,*)'Change ?' +1007 READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 1006 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)'Number of projectiles ? ' +840 READ(5,*,ERR=840)nproj + ELSE + WRITE(*,*)errcom + GOTO 1007 + ENDIF +1006 WRITE(*,'(1x,A11,F7.2)')'Start depth ',x0 + WRITE(*,*)'Change ?' +1009 READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 1008 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)'Start depth ? ' +841 READ(5,'(F7.0)',ERR=841)x0 + ELSE + WRITE(*,*)errcom + GOTO 1009 + ENDIF +1008 WRITE(*,'(1x,A13,F7.2)')'Step size ',cw + WRITE(*,*)'Change ?' +1011 READ(5,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 1010 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,*)' Step size ? ' +842 READ(5,'(F7.0)',ERR=842)cw + ELSE + WRITE(*,*)errcom + GOTO 1011 + ENDIF +1010 continue +c WRITE(*,'(1x,A27,F12.5)')'Random number generator seed',ri +c WRITE(*,*)'Change ?' +c1013 READ(5,'(A1)')f +c IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN +c GOTO 1012 +c ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN +c WRITE(*,*)' Random number seed :' +c843 READ(5,'(F12.0)',ERR=843)ri +c IF((ri/2.-FLOAT(INT(ri/2.))).ne.0.5) ri=ri+1 +c ELSE +c WRITE(*,*)errcom +c GOTO 1013 +c ENDIF +1012 ri2=ri + ri3=ri + WRITE(*,*)' Number of layers ? <7' +850 READ(5,'(I1)',ERR=850)lnum + IF(lnum.LE.0.OR.lnum.GT.7) THEN + WRITE(*,*)errcom + GOTO 1012 + ENDIF +c + CALL NULLEN(ordz,mass,elas,rho,spower,conz, + + conzlayer,thick,ck,sbe,be) + DO l=1,lnum + DO i=1,5 + m=1 +10120 WRITE(*,'(A15,I1)')'Layer number ',l + WRITE(*,'(A15,I1)')'Element number ',m +1014 WRITE(*,'(1x,A44)')'Enter name of element '// + + '(A2), QQ to stop.' +851 READ(5,'(A2)',ERR=851)en + CALL element(en,n) +c + IF(n.EQ.0) THEN + WRITE(*,*)errcom + GOTO 1014 + ENDIF + IF(n.EQ.93) GOTO 1100 + CALL lese(p2,p3,p,n) + ordz(m,l)=FLOAT(n) + mass(m,l)=p2 + elas(m,l)=p3 + sbe(m,l)=30. + DO k=1,5 + spower(k,m,l)=p(k) + ENDDO + WRITE(*,*)' Concentration of an element must be '// + + 'smaller or equal 1.' + WRITE(*,*)' Summe of all concentrations '// + + ' within a layer must be 1.)' +1015 WRITE(*,*)' Enter concentration of the element.' +852 READ(5,'(F7.0)',ERR=852)conz(m,l) + IF(conz(m,l).LE.0.0.OR.conz(m,l).GT.1.0) THEN + WRITE(*,*)errcom + GOTO 1015 + ENDIF + conzlayer(l)=conz(m,l)+conzlayer(l) + IF(conzlayer(l).GT.1.0) THEN + WRITE(*,*)' Sum of concentration within the ' + WRITE(*,'(1x,I1,1x,A29)')l,'. layer is > 1. ' + WRITE(*,*)' S T O P !' + STOP + ENDIF + CALL dichte(p4,n) + rho(l)=rho(l)+p4*conz(m,l) + IF(m.EQ.5) GOTO 1100 +1016 WRITE(*,*)' Another element in layer ?' + READ(5,'(A1)')f + IF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + m=m+1 + GOTO 10120 + ELSEIF(f.EQ.'n'.OR.f.EQ.'N') THEN + ordz(m+1,l)=0.0 + mass(m+1,l)=0.0 + elas(m+1,l)=0.0 + DO k=1,5 + spower(k,m+1,l)=0.0 + ENDDO + IF(m.EQ.1) THEN + CALL dichte(p4,n) + rho(l)=p4 + GOTO 1101 + ENDIF + GOTO 1100 + ELSE + WRITE(*,*)errcom + GOTO 1016 + ENDIF + ENDDO +1100 WRITE(*,*)'Density of layer ? ',rho(l) + WRITE(*,*)'Density OK ?' + READ(5,'(A1)')f + IF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + GOTO 1101 + ELSEIF(f.EQ.'n'.OR.f.EQ.'N') THEN + WRITE(*,*)'Density of layer ? ' +860 READ(5,'(F7.0)',ERR=860)rho(l) + ELSE + WRITE(*,*)errcom + GOTO 1100 + ENDIF +1101 WRITE(*,*)'Thickness of layer ?' +861 READ(5,'(F7.0)',ERR=861)thick(l) + ENDDO + +c WRITE(*,*)' Ausgabe der Layerzusammensetzung auf Bildschirm' +c DO l=1,lnum +c WRITE(*,'(A15,I1)')'Layer Nummer ',l +c DO m=1,5 +c IF(ordz(m,l).EQ.0.0) GOTO 1200 +c WRITE(*,'(A15,I1)')'Element Nummer ',m +c WRITE(*,'(A3,F12.5)')'Z: ',ordz(m,l) +c WRITE(*,'(A3,F12.5)')'M: ',mass(m,l) +c WRITE(*,'(A5,F12.5)')'E-E: ',elas(m,l) +c DO k=1,5 +c WRITE(*,'(A17,I1,1x,F12.5)') +c + 'Stoppkoeffizient ',k,spower(k,m,l) +c ENDDO +c WRITE(*,'(A20,F7.3)')'Konzentration: ',conz(m,l) +c ENDDO +c1200 WRITE(*,'(A43,F7.3)') +c + 'Gesamtmassenkonzentration Layer ->muss 1. sein !',conzlayer(l) +c WRITE(*,'(A30,F7.3)')'Dichte des Layers :',rho(l) +c WRITE(*,'(A30,F12.3)')'Dicke des Layers :',thick(l) +c ENDDO +c +c hier faengt die Schleife fuer die Erstellung verschiedener Dateien +c mit unterschiedlicher Energie und unterschiedlicher Schrittweite an +c + DO 1,n=1,100 +c + IF(n.EQ.1)THEN + IF(batchtype.EQ.'e')THEN + CALL ausgabenam(eproj,parchar) + ELSE +1499 WRITE(*,*)' Specific file name (A3) ' + READ(*,'(A3)',ERR=1499)parchar + ENDIF + ENDIF +c +1500 ausgabe=inpnam//batchtype//parchar//inpext +1501 OPEN(UNIT=11,file=ausgabe,STATUS='NEW',ERR=3000) + CALL batchcreater(inpnam,batchtype,parchar,inpext) +C + WRITE(11,2010)zproj,mproj,eproj,esig,alphaproj,alphaprojsig, + # ef,esb,sheath,erc + WRITE(11,2011)nproj,ri,ri2,ri3,x0,rd,cw,ca,kk0,kk0r,kdee1,kdee2, + # ipot,ipotr,irl + WRITE(11,2012) thick(1),thick(2),thick(3), + # thick(4),thick(5),thick(6),thick(7), + # rho(1),rho(2),rho(3),rho(4),rho(5),rho(6),rho(7), + # ck(1),ck(2),ck(3),ck(4),ck(5),ck(6),ck(7) + DO 2000 I=1,lmax + WRITE(11,2013) ordz(1,I),ordz(2,I),ordz(3,I),ordz(4,I),ordz(5,I) + WRITE(11,2013) mass(1,I),mass(2,I),mass(3,I),mass(4,I),mass(5,I) + WRITE(11,2013) conz(1,I),conz(2,I),conz(3,I),conz(4,I),conz(5,I) + WRITE(11,2013) elas(1,I),elas(2,I),elas(3,I),elas(4,I),elas(5,I) + WRITE(11,2013) sbe(1,I),sbe(2,I),sbe(3,I),sbe(4,I),sbe(5,I) + WRITE(11,2013) be(1,I),be(2,I),be(3,I),be(4,I),be(5,I) +c + WRITE(11,2017) spower(1,1,I),spower(1,2,I),spower(1,3,I), + + spower(1,4,I),spower(1,5,I) + WRITE(11,2017) spower(2,1,I),spower(2,2,I),spower(2,3,I), + + spower(2,4,I),spower(2,5,I) + WRITE(11,2017) spower(3,1,I),spower(3,2,I),spower(3,3,I), + + spower(3,4,I),spower(3,5,I) + WRITE(11,2017) spower(4,1,I),spower(4,2,I),spower(4,3,I), + + spower(4,4,I),spower(4,5,I) + WRITE(11,2017) spower(5,1,I),spower(5,2,I),spower(5,3,I), + + spower(5,4,I),spower(5,5,I) +2000 CONTINUE +c + CLOSE(UNIT=11) +1502 WRITE(*,*)' Another file with new parameters ?' + READ(*,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 9000 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN +c1507 WRITE(*,*)' Welche Parameter sollen geaendert werden ?' +c WRITE(*,*)' E = Energie (default)' +c WRITE(*,*)' S = Esigma' +c WRITE(*,*)' N = Anzahl Myonen' +c WRITE(*,*)' D = Layerdicken' +c READ(*,'(A1)')batchtype + IF(batchtype.EQ.'e'.OR.batchtype.EQ.'E'.or.batchtype.EQ.'') + + THEN + CALL changE(eproj,cw) + CALL ausgabenam(eproj,parchar) + ELSEIF(batchtype.EQ.'s'.OR.batchtype.EQ.'S') THEN + CALL changS(esig) +1520 WRITE(*,*)' spezieller Dateiname (A3) ' + READ(*,'(A3)',ERR=1520)parchar + ELSEIF(batchtype.EQ.'n'.OR.batchtype.EQ.'N') THEN + CALL changN(nproj) +1521 WRITE(*,*)' spezieller Dateiname (A3) ' + READ(*,'(A3)',ERR=1521)parchar + ELSEIF(batchtype.EQ.'d'.OR.batchtype.EQ.'D') THEN + CALL changD(thick,cw,lnum) +1522 WRITE(*,*)' spezieller Dateiname (A3) ' + READ(*,'(A3)',ERR=1522)parchar + ELSEIF(batchtype.EQ.'x'.OR.batchtype.EQ.'X') THEN + CALL changX(eproj,esig) +1523 WRITE(*,*)' spezieller Dateiname (A3) ' + READ(*,'(A3)',ERR=1523)parchar + ELSE + WRITE(*,*)errcom + STOP + ENDIF + GOTO 1500 + ELSE + WRITE(*,*)errcom + GOTO 1502 + ENDIF +3000 WRITE(*,*)' Error during file open. ' + WRITE(*,*)' Enter new file name: ' + READ(*,'(A8)')name + ausgabe=name//inpext + GOTO 1501 +c +1 CONTINUE +c +2010 FORMAT(2F7.2,1F12.2,7F9.2) +2011 FORMAT(I9,3F8.0,1F7.2,1F7.0,2F7.2,6I4,I3) +2012 FORMAT(7F9.2,14F7.2) +2013 FORMAT(5F9.4) +2017 FORMAT(5F12.6) +c +9000 CALL batchexit(inpnam,batchtype) + END +c +c hier werden die Parameter der Elemente eingelesen +c + SUBROUTINE param +c + INTEGER i,n + REAL p2,p3,p4 + REAL p(5) + CHARACTER str2*17,str3*17,str4*69,str5*12 +c Masse der Elemente, Stoppingpowerdaten nach Andersen-Ziegler und Bindungsenergie +c der Elemente im bulk + ENTRY lese(p2,p3,p,n) +c + OPEN(UNIT=31,FILE='masse.dat',STATUS='OLD') + OPEN(UNIT=32,FILE='elast.dat',STATUS='OLD') + OPEN(UNIT=33,FILE='stopicru.dat',STATUS='OLD') +c + DO i=1,n + READ(31,'(A17)')str2 + READ(32,'(A17)')str3 + READ(33,'(A69)')str4 + ENDDO +c + READ(str2(8:16),'(F12.0)')p2 + READ(str3(9:16),'(F12.0)')p3 + READ(str4(4:16),'(F12.0)')p(1) + READ(str4(17:29),'(F12.0)')p(2) + READ(str4(30:42),'(F12.0)')p(3) + READ(str4(43:55),'(F12.0)')p(4) + READ(str4(56:68),'(F12.0)')p(5) +c + CLOSE(UNIT=31) + CLOSE(UNIT=32) + CLOSE(UNIT=33) +c +c OPEN(UNIT=31,POSITION='rewind') +c OPEN(UNIT=32,POSITION='rewind') +c OPEN(UNIT=33,POSITION='rewind') +c + RETURN +c +c die elementaren Dichten +c + ENTRY dichte(p4,n) +c + OPEN(UNIT=34,FILE='dichte.dat',STATUS='OLD') +c + DO i=1,n + READ(34,'(A12)')str5 + ENDDO +c + READ(str5(4:11),'(F12.0)')p4 +c + CLOSE(UNIT=34) +c + RETURN + END +c +c diese Subroutine erkennt das/die Elemente +c + SUBROUTINE element(en,n) +c + INTEGER n + CHARACTER en*2 +c + n=0 +c + IF(en.EQ.' H'.OR.en.EQ.'H') n=1 + IF(en.EQ.'He') n=2 + IF(en.EQ.'Li') n=3 + IF(en.EQ.'Be') n=4 + IF(en.EQ.' B'.OR.en.EQ.'B') n=5 + IF(en.EQ.' C'.OR.en.EQ.'C') n=6 + IF(en.EQ.' N'.OR.en.EQ.'N') n=7 + IF(en.EQ.' O'.OR.en.EQ.'O') n=8 + IF(en.EQ.' F'.OR.en.EQ.'F') n=9 + IF(en.EQ.'Ne') n=10 + IF(en.EQ.'Na') n=11 + IF(en.EQ.'Mg') n=12 + IF(en.EQ.'Al') n=13 + IF(en.EQ.'Si') n=14 + IF(en.EQ.' P'.OR.en.EQ.'P') n=15 + IF(en.EQ.' S'.OR.en.EQ.'S') n=16 + IF(en.EQ.'Cl') n=17 + IF(en.EQ.'Ar') n=18 + IF(en.EQ.' K'.OR.en.EQ.'K') n=19 + IF(en.EQ.'Ca') n=20 + IF(en.EQ.'Sc') n=21 + IF(en.EQ.'Ti') n=22 + IF(en.EQ.' V'.OR.en.EQ.'V') n=23 + IF(en.EQ.'Cr') n=24 + IF(en.EQ.'Mn') n=25 + IF(en.EQ.'Fe') n=26 + IF(en.EQ.'Co') n=27 + IF(en.EQ.'Ni') n=28 + IF(en.EQ.'Cu') n=29 + IF(en.EQ.'Zn') n=30 + IF(en.EQ.'Ga') n=31 + IF(en.EQ.'Ge') n=32 + IF(en.EQ.'As') n=33 + IF(en.EQ.'Se') n=34 + IF(en.EQ.'Br') n=35 + IF(en.EQ.'Kr') n=36 + IF(en.EQ.'Rb') n=37 + IF(en.EQ.'Sr') n=38 + IF(en.EQ.' Y'.OR.en.EQ.'Y') n=39 + IF(en.EQ.'Zr') n=40 + IF(en.EQ.'Nb') n=41 + IF(en.EQ.'Mo') n=42 + IF(en.EQ.'Tc') n=43 + IF(en.EQ.'Ru') n=44 + IF(en.EQ.'Rh') n=45 + IF(en.EQ.'Pd') n=46 + IF(en.EQ.'Ag') n=47 + IF(en.EQ.'Cd') n=48 + IF(en.EQ.'In') n=49 + IF(en.EQ.'Sn') n=50 + IF(en.EQ.'Sb') n=51 + IF(en.EQ.'Te') n=52 + IF(en.EQ.' I'.OR.en.EQ.'I') n=53 + IF(en.EQ.'Xe') n=54 + IF(en.EQ.'Cs') n=55 + IF(en.EQ.'Ba') n=56 + IF(en.EQ.'La') n=57 + IF(en.EQ.'Ce') n=58 + IF(en.EQ.'Pr') n=59 + IF(en.EQ.'Nd') n=60 + IF(en.EQ.'Pm') n=61 + IF(en.EQ.'Sm') n=62 + IF(en.EQ.'Eu') n=63 + IF(en.EQ.'Gd') n=64 + IF(en.EQ.'Tb') n=65 + IF(en.EQ.'Dy') n=66 + IF(en.EQ.'Ho') n=67 + IF(en.EQ.'Er') n=68 + IF(en.EQ.'Tm') n=69 + IF(en.EQ.'Yb') n=70 + IF(en.EQ.'Lu') n=71 + IF(en.EQ.'Hf') n=72 + IF(en.EQ.'Ta') n=73 + IF(en.EQ.' W'.OR.en.EQ.'W') n=74 + IF(en.EQ.'Re') n=75 + IF(en.EQ.'Os') n=76 + IF(en.EQ.'Ir') n=77 + IF(en.EQ.'Pt') n=78 + IF(en.EQ.'Au') n=79 + IF(en.EQ.'Hg') n=80 + IF(en.EQ.'Tl') n=81 + IF(en.EQ.'Pb') n=82 + IF(en.EQ.'Bi') n=83 + IF(en.EQ.'Po') n=84 + IF(en.EQ.'At') n=85 + IF(en.EQ.'Rn') n=86 + IF(en.EQ.'Fr') n=87 + IF(en.EQ.'Ra') n=88 + IF(en.EQ.'Ac') n=89 + IF(en.EQ.'Th') n=90 + IF(en.EQ.'Pa') n=91 + IF(en.EQ.' U'.OR.en.EQ.'U') n=92 + IF(en.EQ.'qq'.OR.en.EQ.'QQ') n=93 +c + RETURN + END +c +c fuers gute Gewissen, die Subroutine nullt alle Parameter +c + SUBROUTINE NULLEN(ordz,mass,elas,rho,spower,conz, + + conzlayer,thick,ck,sbe,be) +c + implicit none + INTEGER i,k,l + REAL*4 ordz(5,7),mass(5,7),elas(5,7),spower(5,5,7),conz(5,7) + REAL*4 conzlayer(7),thick(7),rho(7) + REAL*4 ck(7),sbe(5,7),be(5,7) +c + DO 10 i=1,7 + conzlayer(i)=0.0 + thick(i)=0.0 + rho(i)=0.0 + ck(i)=1.0 + DO 11 k=1,5 + ordz(k,i)=0.0 + mass(k,i)=0.0 + elas(k,i)=0.0 + conz(k,i)=0.0 + sbe(k,i)=0.0 + be(k,i)=0.0 + DO 12 l=1,5 + spower(l,k,i)=0.0 +12 ENDDO +11 ENDDO +10 ENDDO +c + RETURN + END +c +c Subroutine, die die Energie in ein 3zeichigen Character umwandelt +c + SUBROUTINE AUSGABENAM(eproj,parchar) +c + implicit none + + REAL eproj + INTEGER etmp + CHARACTER parchar*3,c1*1,c2*1 +c + DATA c1/'0'/,c2/'0'/ +c + etmp=INT(eproj/1.E2) + IF(etmp.LT.100) THEN + IF(etmp.LT.10) THEN + WRITE(99,'(A1,A1,I1)')c1,c2,etmp + GOTO 10 + ENDIF + WRITE(99,'(A1,I2)')c1,etmp + GOTO 10 + ENDIF + WRITE(99,'(I3)')etmp +10 CLOSE(99) + READ(99,'(A3)')parchar + WRITE(*,'(A3)')parchar + CLOSE(99,STATUS='DELETE') +c + RETURN + END +c +c Subroutine, in der die Parameter Energie, SigmaE, Teilchenanzahl, +c Kanalweite und Schichtdicken geaendert werden koennen +c (nicht aber die Layerzusammensetzung !) +c + SUBROUTINE changes +c + implicit none + + INTEGER i,lnum,nproj + REAL eproj,esig,thick(3),cw + CHARACTER errcom*21,f*1 +c + DATA errcom/'Fehlerhafte Eingabe !'/ +c andere Energie und/oder Kanalweite + ENTRY changE(eproj,cw) + WRITE(*,'(A15,F9.2)')' Old energy: ',eproj +100 WRITE(*,*)' New energy (E < 100.000 eV)' + READ(*,'(F7.0)',ERR=100)eproj +110 WRITE(*,*)' Other step size ?' + READ(*,'(A1)')f + IF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,'(A25,F7.2)')' Old step size: ',cw +115 WRITE(*,*)' New step size' + READ(*,'(F7.0)',ERR=115)cw + ELSEIF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + RETURN + ELSE + WRITE(*,*)errcom + GOTO 110 + ENDIF + RETURN +c andere Energie und anderes Sigma + ENTRY changX(eproj,esig) + WRITE(*,'(A15,F9.2)')' old energy: ',eproj +105 WRITE(*,*)' new energy (E < 100.000 eV)' + READ(*,'(F7.0)',ERR=105)eproj + WRITE(*,'(A15,F7.2)')' old Esigma: ',esig +106 WRITE(*,*)' new Esigma' + READ(*,'(F7.0)',ERR=106)esig + RETURN +c anderes Esigma + ENTRY changS(esig) + WRITE(*,'(A15,F7.2)')' old Esigma: ',esig +120 WRITE(*,*)' new Esigma' + READ(*,'(F7.0)',ERR=120)esig + RETURN +c andere Anzahl Myonen + ENTRY changN(nproj) + WRITE(*,'(A20,I6)')' old number of projectiles: ',nproj +130 WRITE(*,*)' new number of projectiles' + READ(*,'(I6)',ERR=130)nproj + RETURN +c andere Layerdicken und/oder Kanalweite + ENTRY changD(thick,cw,lnum) + WRITE(*,*)lnum + DO i=1,lnum + WRITE(*,'(A23,I1,2x,F7.2)')' old layer thickness ',i, + + thick(i) +140 WRITE(*,'(A23,I1)')' new layer thickness ',i + READ(*,'(F7.0)',ERR=140)thick(i) +150 WRITE(*,*)' Do you want to change the '// + + 'layer thickness ? ' + READ(*,'(A1)')f + IF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + GOTO 160 + ELSEIF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J')THEN + CONTINUE + ELSE + WRITE(*,*)errcom + GOTO 150 + ENDIF + ENDDO +160 WRITE(*,*)' Other step size ?' + READ(*,'(A1)')f + IF(f.EQ.'y'.OR.f.EQ.'Y'.OR.f.EQ.'j'.OR.f.EQ.'J') THEN + WRITE(*,'(A25,F7.2)')' old step size: ',cw + WRITE(*,*)' new step size' + READ(*,'(F7.0)',ERR=130)cw + ELSEIF(f.EQ.'n'.OR.f.EQ.'N'.OR.f.EQ.'') THEN + RETURN + ELSE + WRITE(*,*)errcom + GOTO 160 + ENDIF + RETURN + END +c +c Subroutine, in der die Batchkommandos in die Batch-Datei +c geschrieben werden +c + SUBROUTINE batchcreater(inpnam,batchtype,parchar,inpext) + + implicit none + + CHARACTER inpnam*4,parchar*3,batchtype*1 + CHARACTER inpext*4,rgeext*4,outext*4,batext*4,errext*4 + CHARACTER batchnam*9 +c + DATA rgeext /'.rge'/ + DATA outext /'.out'/ + DATA errext /'.err'/ + DATA batext /'.bat'/ + +100 continue +#if defined (OS_UNIX) + batchnam=inpnam//batchtype +#else + batchnam=inpnam//batchtype//batext +#endif + OPEN(UNIT=19,FILE=batchnam,STATUS='NEW',ERR=900) + GOTO 200 +900 WRITE(*,*)' Batch script already existing, enter new name (A4).' + READ(*,'(A4)')inpnam + GOTO 100 +c +200 continue +#if defined (OS_UNIX) + write(19,'(a)') '#!/bin/csh ' + write(19,'(a)') '#' + WRITE(19,'(a,a,a,a,a,a)') + + 'cp ',inpnam,batchtype,parchar,inpext,' eingabe1.inp' + WRITE(19,'(a)')'trimsp7l_tp' + WRITE(19,'(a,a,a,a,a)') + + 'cp ausgabe1.out ',inpnam,batchtype,parchar,outext + WRITE(19,'(a,a,a,a,a)') + + 'cp ausgabe1.rge ',inpnam,batchtype,parchar,rgeext + WRITE(19,'(a,a,a,a,a)') + + 'cp ausgabe1.err ',inpnam,batchtype,parchar,errext + WRITE(19,'(a)')'rm -f ausgabe1.out' + WRITE(19,'(a)')'rm -f ausgabe1.rge' + WRITE(19,'(a)')'rm -f ausgabe1.err' +#else + WRITE(19,500)'copy',inpnam,batchtype,parchar,inpext, + + 'eingabe1.inp' + WRITE(19,501)'trimsp7l_tp' + WRITE(19,502)'copy ausgabe1.out',inpnam,batchtype,parchar,outext + WRITE(19,502)'copy ausgabe1.rge',inpnam,batchtype,parchar,rgeext + WRITE(19,502)'copy ausgabe1.err',inpnam,batchtype,parchar,errext + WRITE(19,505)'del ausgabe1.out' + WRITE(19,505)'del ausgabe1.rge' + WRITE(19,505)'del ausgabe1.err' +#endif + RETURN +c + ENTRY batchexit(inpnam,batchtype) +#if defined (OS_UNIX) + WRITE(19,'(a,a,a,a)')'mv fort.33 ',inpnam,batchtype,'.dat' + WRITE(19,'(a)')'rm -f eingabe1.inp' + WRITE(19,'(a)')'rm -f edist' +#else + WRITE(19,510)'ren fort.33',inpnam,batchtype,'.dat' + WRITE(19,505)'del eingabe1.inp' + WRITE(19,506)'del edist' + WRITE(19,511)'exit' +#endif + RETURN +c +500 FORMAT(A4,1x,A4,A1,A3,A4,1x,A12) +501 FORMAT(A8) +502 FORMAT(A17,1x,A4,A1,A3,A4) +505 FORMAT(A16) +506 FORMAT(A9) +510 FORMAT(A11,1x,A4,A1,A4) +511 FORMAT(A4) + END + + + diff --git a/trimsp/obsolete/trimsp7l_tp.F b/trimsp/obsolete/trimsp7l_tp.F index afccbd7..2ed4af5 100644 --- a/trimsp/obsolete/trimsp7l_tp.F +++ b/trimsp/obsolete/trimsp7l_tp.F @@ -5225,3 +5225,298 @@ C RETURN END + SUBROUTINE RANLUX(RVEC,LENV) +C Subtract-and-borrow random number generator proposed by +C Marsaglia and Zaman, implemented by F. James with the name +C RCARRY in 1991, and later improved by Martin Luescher +C in 1993 to produce "Luxury Pseudorandom Numbers". +C Fortran 77 coded by F. James, 1993 +C +C LUXURY LEVELS. +C ------ ------ The available luxury levels are: +C +C level 0 (p=24): equivalent to the original RCARRY of Marsaglia +C and Zaman, very long period, but fails many tests. +C level 1 (p=48): considerable improvement in quality over level 0, +C now passes the gap test, but still fails spectral test. +C level 2 (p=97): passes all known tests, but theoretically still +C defective. +C level 3 (p=223): DEFAULT VALUE. Any theoretically possible +C correlations have very small chance of being observed. +C level 4 (p=389): highest possible luxury, all 24 bits chaotic. +C +C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +C!!! Calling sequences for RANLUX: ++ +C!!! CALL RANLUX (RVEC, LEN) returns a vector RVEC of LEN ++ +C!!! 32-bit random floating point numbers between ++ +C!!! zero (not included) and one (also not incl.). ++ +C!!! CALL RLUXGO(LUX,INT,K1,K2) initializes the generator from ++ +C!!! one 32-bit integer INT and sets Luxury Level LUX ++ +C!!! which is integer between zero and MAXLEV, or if ++ +C!!! LUX .GT. 24, it sets p=LUX directly. K1 and K2 ++ +C!!! should be set to zero unless restarting at a break++ +C!!! point given by output of RLUXAT (see RLUXAT). ++ +C!!! CALL RLUXAT(LUX,INT,K1,K2) gets the values of four integers++ +C!!! which can be used to restart the RANLUX generator ++ +C!!! at the current point by calling RLUXGO. K1 and K2++ +C!!! specify how many numbers were generated since the ++ +C!!! initialization with LUX and INT. The restarting ++ +C!!! skips over K1+K2*E9 numbers, so it can be long.++ +C!!! A more efficient but less convenient way of restarting is by: ++ +C!!! CALL RLUXIN(ISVEC) restarts the generator from vector ++ +C!!! ISVEC of 25 32-bit integers (see RLUXUT) ++ +C!!! CALL RLUXUT(ISVEC) outputs the current values of the 25 ++ +C!!! 32-bit integer seeds, to be used for restarting ++ +C!!! ISVEC must be dimensioned 25 in the calling program ++ +C!!! ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + DIMENSION RVEC(LENV) + DIMENSION SEEDS(24), ISEEDS(24), ISDEXT(25) + PARAMETER (MAXLEV=4, LXDFLT=3) + DIMENSION NDSKIP(0:MAXLEV) + DIMENSION NEXT(24) + PARAMETER (TWOP12=4096., IGIGA=1000000000,JSDFLT=314159265) + PARAMETER (ITWO24=2**24, ICONS=2147483563) + SAVE NOTYET, I24, J24, CARRY, SEEDS, TWOM24, TWOM12, LUXLEV + SAVE NSKIP, NDSKIP, IN24, NEXT, KOUNT, MKOUNT, INSEED + INTEGER LUXLEV + LOGICAL NOTYET + DATA NOTYET, LUXLEV, IN24, KOUNT, MKOUNT /.TRUE., LXDFLT, 0,0,0/ + DATA I24,J24,CARRY/24,10,0./ +C default +C Luxury Level 0 1 2 *3* 4 + DATA NDSKIP/0, 24, 73, 199, 365 / +Corresponds to p=24 48 97 223 389 +C time factor 1 2 3 6 10 on slow workstation +C 1 1.5 2 3 5 on fast mainframe +C +C NOTYET is .TRUE. if no initialization has been performed yet. +C Default Initialization by Multiplicative Congruential + IF (NOTYET) THEN + NOTYET = .FALSE. + JSEED = JSDFLT + INSEED = JSEED + LUXLEV = LXDFLT + NSKIP = NDSKIP(LUXLEV) + LP = NSKIP + 24 + IN24 = 0 + KOUNT = 0 + MKOUNT = 0 + TWOM24 = 1. + DO 25 I= 1, 24 + TWOM24 = TWOM24 * 0.5 + K = JSEED/53668 + JSEED = 40014*(JSEED-K*53668) -K*12211 + IF (JSEED .LT. 0) JSEED = JSEED+ICONS + ISEEDS(I) = MOD(JSEED,ITWO24) + 25 CONTINUE + TWOM12 = TWOM24 * 4096. + DO 50 I= 1,24 + SEEDS(I) = REAL(ISEEDS(I))*TWOM24 + NEXT(I) = I-1 + 50 CONTINUE + NEXT(1) = 24 + I24 = 24 + J24 = 10 + CARRY = 0. + IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 + ENDIF +C +C The Generator proper: "Subtract-with-borrow", +C as proposed by Marsaglia and Zaman, +C Florida State University, March, 1989 +C + DO 100 IVEC= 1, LENV + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + RVEC(IVEC) = UNI +C small numbers (with less than 12 "significant" bits) are "padded". + IF (UNI .LT. TWOM12) THEN + RVEC(IVEC) = RVEC(IVEC) + TWOM24*SEEDS(J24) +C and zero is forbidden in case someone takes a logarithm + IF (RVEC(IVEC) .EQ. 0.) RVEC(IVEC) = TWOM24*TWOM24 + ENDIF +C Skipping to luxury. As proposed by Martin Luscher. + IN24 = IN24 + 1 + IF (IN24 .EQ. 24) THEN + IN24 = 0 + KOUNT = KOUNT + NSKIP + DO 90 ISK= 1, NSKIP + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + 90 CONTINUE + ENDIF + 100 CONTINUE + KOUNT = KOUNT + LENV + IF (KOUNT .GE. IGIGA) THEN + MKOUNT = MKOUNT + 1 + KOUNT = KOUNT - IGIGA + ENDIF + RETURN +C +C Entry to input and float integer seeds from previous run + ENTRY RLUXIN(ISDEXT) + NOTYET = .FALSE. + TWOM24 = 1. + DO 195 I= 1, 24 + NEXT(I) = I-1 + 195 TWOM24 = TWOM24 * 0.5 + NEXT(1) = 24 + TWOM12 = TWOM24 * 4096. + WRITE(6,'(A)') ' FULL INITIALIZATION OF RANLUX WITH 25 INTEGERS:' + WRITE(6,'(5X,5I12)') ISDEXT + DO 200 I= 1, 24 + SEEDS(I) = REAL(ISDEXT(I))*TWOM24 + 200 CONTINUE + CARRY = 0. + IF (ISDEXT(25) .LT. 0) CARRY = TWOM24 + ISD = IABS(ISDEXT(25)) + I24 = MOD(ISD,100) + ISD = ISD/100 + J24 = MOD(ISD,100) + ISD = ISD/100 + IN24 = MOD(ISD,100) + ISD = ISD/100 + LUXLEV = ISD + IF (LUXLEV .LE. MAXLEV) THEN + NSKIP = NDSKIP(LUXLEV) + WRITE (6,'(A,I2)') ' RANLUX LUXURY LEVEL SET BY RLUXIN TO: ', + + LUXLEV + ELSE IF (LUXLEV .GE. 24) THEN + NSKIP = LUXLEV - 24 + WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXIN TO:',LUXLEV + ELSE + NSKIP = NDSKIP(MAXLEV) + WRITE (6,'(A,I5)') ' RANLUX ILLEGAL LUXURY RLUXIN: ',LUXLEV + LUXLEV = MAXLEV + ENDIF + INSEED = -1 + RETURN +C +C Entry to ouput seeds as integers + ENTRY RLUXUT(ISDEXT) + DO 300 I= 1, 24 + ISDEXT(I) = INT(SEEDS(I)*TWOP12*TWOP12) + 300 CONTINUE + ISDEXT(25) = I24 + 100*J24 + 10000*IN24 + 1000000*LUXLEV + IF (CARRY .GT. 0.) ISDEXT(25) = -ISDEXT(25) + RETURN +C +C Entry to output the "convenient" restart point + ENTRY RLUXAT(LOUT,INOUT,K1,K2) + LOUT = LUXLEV + INOUT = INSEED + K1 = KOUNT + K2 = MKOUNT + RETURN +C +C Entry to initialize from one or three integers + ENTRY RLUXGO(LUX,INS,K1,K2) + IF (LUX .LT. 0) THEN + LUXLEV = LXDFLT + ELSE IF (LUX .LE. MAXLEV) THEN + LUXLEV = LUX + ELSE IF (LUX .LT. 24 .OR. LUX .GT. 2000) THEN + LUXLEV = MAXLEV + WRITE (6,'(A,I7)') ' RANLUX ILLEGAL LUXURY RLUXGO: ',LUX + ELSE + LUXLEV = LUX + DO 310 ILX= 0, MAXLEV + IF (LUX .EQ. NDSKIP(ILX)+24) LUXLEV = ILX + 310 CONTINUE + ENDIF + IF (LUXLEV .LE. MAXLEV) THEN + NSKIP = NDSKIP(LUXLEV) + WRITE(6,'(A,I2,A,I4)') ' RANLUX LUXURY LEVEL SET BY RLUXGO :', + + LUXLEV,' P=', NSKIP+24 + ELSE + NSKIP = LUXLEV - 24 + WRITE (6,'(A,I5)') ' RANLUX P-VALUE SET BY RLUXGO TO:',LUXLEV + ENDIF + IN24 = 0 + IF (INS .LT. 0) WRITE (6,'(A)') + + ' Illegal initialization by RLUXGO, negative input seed' + IF (INS .GT. 0) THEN + JSEED = INS + WRITE(6,'(A,3I12)') ' RANLUX INITIALIZED BY RLUXGO FROM SEEDS', + + JSEED, K1,K2 + ELSE + JSEED = JSDFLT + WRITE(6,'(A)')' RANLUX INITIALIZED BY RLUXGO FROM DEFAULT SEED' + ENDIF + INSEED = JSEED + NOTYET = .FALSE. + TWOM24 = 1. + DO 325 I= 1, 24 + TWOM24 = TWOM24 * 0.5 + K = JSEED/53668 + JSEED = 40014*(JSEED-K*53668) -K*12211 + IF (JSEED .LT. 0) JSEED = JSEED+ICONS + ISEEDS(I) = MOD(JSEED,ITWO24) + 325 CONTINUE + TWOM12 = TWOM24 * 4096. + DO 350 I= 1,24 + SEEDS(I) = REAL(ISEEDS(I))*TWOM24 + NEXT(I) = I-1 + 350 CONTINUE + NEXT(1) = 24 + I24 = 24 + J24 = 10 + CARRY = 0. + IF (SEEDS(24) .EQ. 0.) CARRY = TWOM24 +C If restarting at a break point, skip K1 + IGIGA*K2 +C Note that this is the number of numbers delivered to +C the user PLUS the number skipped (if luxury .GT. 0). + KOUNT = K1 + MKOUNT = K2 + IF (K1+K2 .NE. 0) THEN + DO 500 IOUTER= 1, K2+1 + INNER = IGIGA + IF (IOUTER .EQ. K2+1) INNER = K1 + DO 450 ISK= 1, INNER + UNI = SEEDS(J24) - SEEDS(I24) - CARRY + IF (UNI .LT. 0.) THEN + UNI = UNI + 1.0 + CARRY = TWOM24 + ELSE + CARRY = 0. + ENDIF + SEEDS(I24) = UNI + I24 = NEXT(I24) + J24 = NEXT(J24) + 450 CONTINUE + 500 CONTINUE +C Get the right value of IN24 by direct calculation + IN24 = MOD(KOUNT, NSKIP+24) + IF (MKOUNT .GT. 0) THEN + IZIP = MOD(IGIGA, NSKIP+24) + IZIP2 = MKOUNT*IZIP + IN24 + IN24 = MOD(IZIP2, NSKIP+24) + ENDIF +C Now IN24 had better be between zero and 23 inclusive + IF (IN24 .GT. 23) THEN + WRITE (6,'(A/A,3I11,A,I5)') + + ' Error in RESTARTING with RLUXGO:',' The values', INS, + + K1, K2, ' cannot occur at luxury level', LUXLEV + IN24 = 0 + ENDIF + ENDIF + RETURN + END + +