programma di Enrico recuperato da Shimizu

ciao

silvia

-------------------------- Messaggio originale ---------------------------
Oggetto: GOE program
Da:      "Yoshifumi R. Shimizu" <yrsh2scp@mbox.nc.kyushu-u.ac.jp>
Data:    Mar, 26 Settembre 2006 4:15 am
A:       silvia.leoni@mi.infn.it
Cc:      matsuo@nt.sc.niigata-u.ac.jp
         yrsh2scp@mbox.nc.kyushu-u.ac.jp
--------------------------------------------------------------------------

Dear Silvia:

     This is Shimizu, so long.  It looks that you and people in Milano
are working hard, and getting very interesting data on excited SD bands.
I'm very happy to contribute.

     About the inquiry of yours,

>At the moment, to calculate the N_out factor we use the approximate
expression given in ref. Shimizu et al., NPA557 (1993)99c (eq. 3.4).
Discussing with Enrico we realized that this is correct only for small
values of N_out (tipically < 1/2).

Yes.  Enrico sent me the program long time ago, and we are using it. I
don't know which form is convenient for you, so I just send the
original one (I think quite near to the original, I'm not so sure). But I
must say, Enrico should have found it on his computer!

      I just briefly explain the meaning of inputs, though Enrico can
tell you about them.

---------------------------------------------------------------
        read*,nstep,nn
        read*,gdLOG1,gdLOG2,gdint
        read*,g2g1LOG1,g2g1LOG2,g2g1int
        READ*,INAG,ISTORE,IFORMAT,ICONTR,ICONTR1,icontrv
---------------------------------------------------------------
C 1000,40      /nstep,nn$
C -4,1,0.1     /gdLOG1,gdLOG2,gdint$
C -2,2,0.1     /g2g1LOG1,g2g1LOG2,g2g1int$
C 0 0 0 1 0 0  /INAG,ISTORE,IFORMAT,ICONTR,ICONTR1,icontrv$
---------------------------------------------------------------

nstep     --- number of simulations for average
nn        --- dimension of GOE matrix
gdLOG1    --- initial Gamma/D
gdLOG2    --- final Gamma/D
gdint     --- mesh interval of Gamma/D
g2g1LOG1  --- initial Gamma_SD/Gamma_ND
g2g1LOG2  --- final Gamma_SD/Gamma_ND
g2g1int   --- mesh interval of Gamma_SD/Gamma_ND

For other control parameters, INAG etc, I don't know their precise
meanings; just use the values of the example input.

The output is written on units 13, and 14.  Those in 13 is F_out, and 14
is it's variance.  So for usual calculations output of only 13 is
necessary. The output are large numbers of data according to the input
mesh of Gamma/D and Gamma_SD/Gamma_ND.

Oh, the most important is that all the quantities, Gamma/D,
Gamma_SD/Gamma_ND, and F_out, and variance are in LOG10; so you have to
make 10^(x) before using them.

     Actually the program is not so difficult to see what is what.  But,
if you have any questions, don't hesitate to ask me.

     Please send best regards to people in Milano, Angera, Ricardo, Enrico,
Pia-Franchesco ...

     Best regards,              Yoshifumi Shimizu
------------------------------------------------------------------------
input file
------------------------------------------------------------------------
1000  , 40     /Nstep,Ndim
-4.0 2.0  0.1  /GDlog1,GDlog2,GDint
-2.0 2.0  0.1  /G2G1log1,G2G1log2,G2G1int
0 0 0 1 0 0    /Inag,Istore,Iformat,Icontr,Icontr1,Icontrv
------------------------------------------------------------------------
fortran source program
------------------------------------------------------------------------ C
$ ASS P.OUT FOR013
C $ R SH
C 1000,40      /nstep,nn
C -4,1,0.1     /gdLOG1,gdLOG2,gdint
C -2,2,0.1     /g2g1LOG1,g2g1LOG2,g2g1int
C 0 0 0 1 0 0  /INAG,ISTORE,IFORMAT,ICONTR,ICONTR1,icontrv
C
        PARAMETER(NTRY=1000,NDIM=1000,NDATA=100)
        COMMON  VV(NDIM,NDIM),EE(NDIM),DD(NDIM)
        common  streng(NTRY),pp(NTRY,NDIM)
        common  pmat(NDATA),VARMAT(NDATA)
        character*50 file_name
        COMMON/CONTR/ICONTR,ICONTR1,icontrv
        PI=4.*ATAN(1.)
        read*,nstep,nn
        read*,gdLOG1,gdLOG2,gdint
        read*,g2g1LOG1,g2g1LOG2,g2g1int
        READ*,INAG,ISTORE,IFORMAT,ICONTR,ICONTR1,icontrv
        write(6,*) ' nstep = ',nstep, ' nn = ',nn

c
        if(nn.gt.NDIM) then
        write(6,*)'Input dim (nn) exceeds the limit (NMAX), stop!' stop
        endif
        if(nstep.gt.NTRY) then
        write(6,*)'Input dim (nstep) exceeds the limit (NTRY), stop!' stop
        endif
c

        WRITE(13,1003) gdlog1,gdlog2,GDINT
        WRITE(13,1003) g2g1LOG1,g2g1LOG2,g2g1INT
        WRITE(14,1003) gdlog1,gdlog2,GDINT
        WRITE(14,1003) g2g1LOG1,g2g1LOG2,g2g1INT
1003  FORMAT(2X,3(F9.5,1X))
c
        i=(g2g1LOG2-g2g1LOG1)/g2g1INT+1.5
        if(i.gt.NDATA) then
        write(6,*)'Input mesh data (g2g1) exceeds the limit, stop!' stop
        endif

        do 18 glog= gdlog1,gdlog2+1.e-6,gdint

        gdrefer= 10.**(glog)
        gd=gdrefer
        dist=1.
        gamma=gd
        do 10 in=1,nstep
c        write(6,*) ' glog,i, gamma '
c        write(6,*) glog,i,gamma
c         if(inag.ne.0) eps1= x02ade(it)
        SIGV=SQRT(GAMMA*DIST/(2.*PI))
        CALL SETMAT(NN,VV,SIGV,DIST)
         IF(INAG.EQ.0)CALL TRED2(VV,NN,DD,EE)
         IF(INAG.EQ.0)CALL TQLI(DD,EE,NN,VV)
c        IF(INAG.NE.0) CALL F01AJE(NN,EPS1,VV,NDIM,D,E,V,NDIM)
        IFAIL=0
c        IF(INAG.NE.0) CALL F02AME(NN,EPS2,D,E,V,NDIM,IFAIL)
        IF(IFAIL.NE.0) WRITE(6,*) IFAIL
        IF(IFAIL.NE.0) STOP
        DO 300 J=1,NN
        IF(INAG.EQ.0)STRENG(J)=VV(1,J)**2

C       IF(INAG.NE.0)STRENG(J)=V(1,J)**2
        PP(IN,J) = STRENG(J)

 300    continue
        IF(IFORMAT.EQ.0.AND.ISTORE.EQ.1)
     1  write(12) (STRENG(J),J=1,NN)
        IF(IFORMAT.NE.0.AND.ISTORE.EQ.1)
     1  write(12,*)(STRENG(J),J=1,NN)
 10     continue

        WRITE(13,1002)gLOG, GD
        WRITE(14,1002)gLOG, GD
1002  FORMAT(5X,F9.5,1X,E12.5)
        ICOUNT=1
        DO 90 g2g1LOG=g2g1LOG1,g2g1LOG2+1.e-6,g2g1int
        ag2g1refer= 10.**(g2g1log)
        ag2g1=ag2g1refer
        pout=0.
        SQOUT=0.
        do 11 in= 1,nstep
        PLEAVE=0.
        DO 310 J=1,NN
        STREN=PP(IN,J)
        S1=(1-STREN)
        S2=STREN
        SFAC= S1/(S1+S2*Ag2g1)
        PLEAVE=PLEAVE+ STREN*SFAC
 310    continue
        POUT= PLEAVE +POUT
        SQOUT= PLEAVE*PLEAVE+SQOUT
 11     CONTINUE
        POUT=POUT/NSTEP
        SQOUT=SQOUT/NSTEP
        PMAT(ICOUNT)=POUT
        VARMAT(ICOUNT)= SQOUT- POUT*POUT
        ICOUNT=ICOUNT+1
 90     CONTINUE
c        write(13,1001) (pmat(i),I=1,ICOUNT-1)
c        write(14,1001) (VARmat(i),i=1,ICOUNT-1)
        write(13,1001) (log10(pmat(i)),I=1,ICOUNT-1)
        write(14,1001) (log10(VARmat(i)),i=1,ICOUNT-1)
1001  FORMAT(5(F9.5,1X))
        if(istore.eq.1)close(12)

 18     continue


 510       END

C **********************************************************************

        SUBROUTINE SETMAT(NN,VV,SIGV,DIST)
        PARAMETER(NDIM=1000)
        DOUBLE PRECISION DSEED
        DATA DSEED/123457.D0/
        COMMON/CONTR/ICONTR,ICONTR1,icontrv
        DIMENSION VV(NDIM,*)
        PI=4.*ATAN(1.)
C ...    S. D. STATE AT ENERGY 0
        VV(1,1)=0.
C ...    S. D. COUPLING
        DO 100 I=2,NN
        CHI1=RANDOM(DSEED)
        CHI2=RANDOM(DSEED)
        if(icontrv.eq.0) then
        IF(ICONTR.EQ.0) THEN
        U=ALOG(1-CHI1)
        VV(1,I)=COS(CHI2*PI)*SIGV*U
        ELSE
        U=-2.*ALOG(1-CHI1)
        VV(1,I)=COS(2.*CHI2*PI)*SIGV*SQRT(U)
        END IF
        else
        vv(1,I)= sigv*sqrt(2./pi)
        endif
        VV(I,1)=VV(1,I)
  100   CONTINUE
C ...    G. O. E. COUPLING
        IF(ICONTR1.EQ.0)VGOE=SQRT(NN-1.)*DIST/PI
        IF(ICONTR1.NE.0)VGOE=SQRT(NN-1.)*DIST/(0.9*SQRT(2.)*PI)
        NN1=NN-1
        DO 200 I=2,NN1
        DO 200 J=I,NN
        CHI1=RANDOM(DSEED)
        CHI2=RANDOM(DSEED)
        IF(ICONTR.EQ.0) THEN
        U=ALOG(1-CHI1)
        VV(I,J)=COS(CHI2*PI)*VGOE*U
        ELSE
        U=-2.*ALOG(1-CHI1)
        VV(I,J)=COS(2.*CHI2*PI)*VGOE*SQRT(U)
        ENDIF
  200   VV(J,I)=VV(I,J)
        DO 250 I=2,NN
        CHI1=RANDOM(DSEED)
        CHI2=RANDOM(DSEED)
        IF(ICONTR.EQ.0) THEN
        U=ALOG(1-CHI1)
        VV(I,I)=SQRT(2.)*COS(CHI2*PI)*VGOE*U
        ELSE
        U=-2.*ALOG(1-CHI1)
        VV(I,I)=SQRT(2.)*COS(2.*CHI2*PI)*VGOE*SQRT(U)
        ENDIF
 250    CONTINUE
        RETURN
        END

C **********************************************************************

        FUNCTION RANDOM(DSEED)
        DOUBLE PRECISION DSEED
        DOUBLE PRECISION D2P31M,D2P31
        DATA D2P31M /2147483647.D0 /
        DATA D2P31  /2147483711.D0 /
        DSEED=DMOD(16807.D0*DSEED,D2P31M)
        RANDOM=DSEED/D2P31
        RETURN
        END

C **********************************************************************

      SUBROUTINE TRED2(A,NV,D,E)
      PARAMETER(NDIM=1000)
      DIMENSION A(NDIM,*),D(*),E(*)

      IF(NV.GT.1)THEN
        DO 18 I=NV,2,-1
          L=I-1
          H=0.
          SCALE=0.
          IF(L.GT.1)THEN
            DO 11 K=1,L
              SCALE=SCALE+ABS(A(I,K))
11          CONTINUE
            IF(SCALE.EQ.0.)THEN
              E(I)=A(I,L)
            ELSE
              DO 12 K=1,L
                A(I,K)=A(I,K)/SCALE
                H=H+A(I,K)**2
12            CONTINUE
              F=A(I,L)
              G=-SIGN(SQRT(H),F)
              E(I)=SCALE*G
              H=H-F*G
              A(I,L)=F-G
              F=0.
              DO 15 J=1,L
                A(J,I)=A(I,J)/H
                G=0.
                DO 13 K=1,J
                  G=G+A(J,K)*A(I,K)
13              CONTINUE
                IF(L.GT.J)THEN
                  DO 14 K=J+1,L
                    G=G+A(K,J)*A(I,K)
14                CONTINUE
                ENDIF
                E(J)=G/H
                F=F+E(J)*A(I,J)
15            CONTINUE
              HH=F/(H+H)
              DO 17 J=1,L
                F=A(I,J)
                G=E(J)-HH*F
                E(J)=G
                DO 16 K=1,J
                  A(J,K)=A(J,K)-F*E(K)-G*A(I,K)
16              CONTINUE
17            CONTINUE
            ENDIF
          ELSE
            E(I)=A(I,L)
          ENDIF
          D(I)=H
18      CONTINUE
      ENDIF
      D(1)=0.
      E(1)=0.
      DO 23 I=1,NV
        L=I-1
        IF(D(I).NE.0.)THEN
          DO 21 J=1,L
            G=0.
            DO 19 K=1,L
              G=G+A(I,K)*A(K,J)
19          CONTINUE
            DO 20 K=1,L
              A(K,J)=A(K,J)-G*A(K,I)
20          CONTINUE
21        CONTINUE
        ENDIF
        D(I)=A(I,I)
        A(I,I)=1.
        IF(L.GE.1)THEN
          DO 22 J=1,L
            A(I,J)=0.
            A(J,I)=0.
22        CONTINUE
        ENDIF
23    CONTINUE

      RETURN
      END

C **********************************************************************

      SUBROUTINE TQLI(D,E,NV,Z)
      PARAMETER(NDIM=1000)
      DIMENSION D(*),E(*),Z(NDIM,*)
      IF (NV.GT.1) THEN
        DO 11 I=2,NV
          E(I-1)=E(I)
11      CONTINUE
        E(NV)=0.
        DO 15 L=1,NV
          ITER=0
1         DO 12 M=L,NV-1
            DD=ABS(D(M))+ABS(D(M+1))
            IF (ABS(E(M))+DD.EQ.DD) GO TO 2
12        CONTINUE
          M=NV
2         IF(M.NE.L)THEN
            IF(ITER.EQ.30)PAUSE 'too many iterations'
            ITER=ITER+1
            G=(D(L+1)-D(L))/(2.*E(L))
            R=SQRT(G**2+1.)
            G=D(M)-D(L)+E(L)/(G+SIGN(R,G))
            S=1.
            C=1.
            P=0.
            DO 14 I=M-1,L,-1
              F=S*E(I)
              B=C*E(I)
              IF(ABS(F).GE.ABS(G))THEN
                C=G/F
                R=SQRT(C**2+1.)
                E(I+1)=F*R
                S=1./R
                C=C*S
              ELSE
                S=F/G
                R=SQRT(S**2+1.)
                E(I+1)=G*R
                C=1./R
                S=S*C
              ENDIF
              G=D(I+1)-P
              R=(D(I)-G)*S+2.*C*B
              P=S*R
              D(I+1)=G+P
              G=C*R-B
              DO 13 K=1,NV
                F=Z(K,I+1)
                Z(K,I+1)=S*Z(K,I)+C*F
                Z(K,I)=C*Z(K,I)-S*F
13            CONTINUE
14          CONTINUE
            D(L)=D(L)-P
            E(L)=G
            E(M)=0.
            GO TO 1
          ENDIF
15      CONTINUE
      ENDIF
      RETURN
      END



-- 
Silvia Leoni
Dipartimento di Fisica
Via Celoria 16
20133 Milano, Italy
Tel. +39 0250317261
Fax. +39 0250317487
web-page: http://www.mi.infn.it/~sleoni
