**Listed On  Mon Jan 14 18:56:35 1991
  **
**List of Fold**       "proptest"
**List of File**       "proptest.tsr"
**File Last Modified    Wed Nov  7 09:43:03 1990
**List all lines  with Fold Headers 
**Excluding : NO LIST folds 

-- Exact calculation of some free-fermion propagators with
-- smeared operators.

#USE biolib

CHAN OF ANY to.save :
PAR
  screen IS to.save :
  SEQ
    --{{{  
    
    VAL l.t IS 8 :
    VAL l.s IS 4 :
    
    --{{{  PROC fermion.propagator ([] REAL32 ^ 5, VAL REAL32 ^ 2, VAL INT ^ 2)
    PROC fermion.propagator ([] [] [] [] REAL32 s.1, s.t, s.x, s.y, s.z,
                            VAL REAL32 k.fermion, r.fermion, VAL INT r.ini, r.fin)
    
      -- Fermion propagator, by gamma matrix components.
    
      -- Result is given as mixed real-space/momentum-space function -
      -- time dependence in real-space, spatial dependence in momentum-
      -- space.
    
      -- 's.1', 's.t' are pure real.  The rest are pure imaginary.
    
      --{{{  trig tables
      
      [l.t] [l.t] REAL32 re.phase, im.phase :
      [l.t] REAL32 cos.t, sin.t :
      [l.s] REAL32 cos.s, sin.s :
      [l.s] REAL32 v.s :
      
      --}}}  
    
      --{{{  PROC trig.tables ()
      PROC trig.tables ()
      
        -- Initialise trig tables (and related stuff).
      
        #USE snglmath
      
        VAL pi     IS 3.1415926 (REAL32) :
        VAL unit.t IS (pi + pi) / (REAL32 ROUND l.t) :
        VAL unit.s IS (pi + pi) / (REAL32 ROUND l.s) :
      
        --{{{  REAL32 FUNCTION v (VAL INT, VAL REAL32)
        REAL32 FUNCTION v (VAL INT r, VAL REAL32 p)
        
          -- "Smearing" kernel function in momentum space.
          -- 'r' is smearing radius.
        
          REAL32 v :
          VALOF
            IF
              p = 0.0 (REAL32)
                v := 1.0 (REAL32)
              TRUE
                VAL d      IS REAL32 ROUND ((2 * r) + 1) :
                VAL p.by.2 IS p / 2.0 (REAL32) :
                v := SIN (d * p.by.2) / (d * SIN (p.by.2))
            RESULT v
        :
        --}}}  
      
        SEQ
          --{{{  "phase" for time Fourier transform
          SEQ t = 0 FOR l.t
            SEQ k = 0 FOR l.t
              VAL arg IS (REAL32 ROUND ((k * t) REM l.t)) * unit.t :
              SEQ
                re.phase [t] [k] := COS (arg)
                im.phase [t] [k] := SIN (arg)
          --}}}  
          --{{{  Trig functions of k
          SEQ
            SEQ k = 0 FOR l.t
              VAL arg IS (REAL32 ROUND k) * unit.t :
              SEQ
                cos.t [k] := COS (arg)
                sin.t [k] := SIN (arg)
            SEQ k = 0 FOR l.s
              VAL arg IS (REAL32 ROUND k) * unit.s :
              SEQ
                cos.s [k] := COS (arg)
                sin.s [k] := SIN (arg)
          --}}}  
          --{{{  Smearing kernel
          SEQ k = 0 FOR l.s
            VAL p IS (REAL32 ROUND k) * unit.s :
            v.s [k] := v (r.ini, p) * v (r.fin, p)
          --}}}  
      :
      --}}}  
    
      [l.t] [l.s] [l.s] [l.s] REAL32 r.1, r.t, r.x, r.y, r.z :
      SEQ
        trig.tables ()
        --{{{  First calculate 'r' - pure momentum-space version
        
        -- 'r.1' is real.  Rest are imaginary.
        
        SEQ k.t = 0 FOR l.t
          VAL v IS 1.0 (REAL32) :
          SEQ k.x = 0 FOR l.s
            VAL v IS v.s [k.x] * v :
            SEQ k.y = 0 FOR l.s
              VAL v IS v.s [k.y] * v :
              SEQ k.z = 0 FOR l.s
                VAL v IS v.s [k.z] * v :
                --{{{  VAL sin's and cos's
                
                VAL sin.k.t IS sin.t [k.t] :
                VAL sin.k.x IS sin.s [k.x] :
                VAL sin.k.y IS sin.s [k.y] :
                VAL sin.k.z IS sin.s [k.z] :
                
                VAL cos.k.t IS cos.t [k.t] :
                VAL cos.k.x IS cos.s [k.x] :
                VAL cos.k.y IS cos.s [k.y] :
                VAL cos.k.z IS cos.s [k.z] :
                
                --}}}  
                --{{{  VAL numerator terms
                
                VAL two.k  IS 2.0 (REAL32) * k.fermion :
                VAL two.kr IS two.k * r.fermion :
                
                VAL num.1  IS 1.0 (REAL32) -
                  (two.kr * (cos.k.t + (cos.k.x + (cos.k.y + cos.k.z)))) :
                
                VAL num.t  IS two.k * sin.k.t :
                VAL num.x  IS two.k * sin.k.x :
                VAL num.y  IS two.k * sin.k.y :
                VAL num.z  IS two.k * sin.k.z :
                
                --}}}  
                VAL den IS (num.1 * num.1) + ((num.t * num.t) +
                          ((num.x * num.x) + ((num.y * num.y) + (num.z * num.z)))) :
                SEQ
                  r.1 [k.t] [k.x] [k.y] [k.z] := v * (num.1 / den)
                  r.t [k.t] [k.x] [k.y] [k.z] := v * (num.t / den)
                  r.x [k.t] [k.x] [k.y] [k.z] := v * (num.x / den)
                  r.y [k.t] [k.x] [k.y] [k.z] := v * (num.y / den)
                  r.z [k.t] [k.x] [k.y] [k.z] := v * (num.z / den)
        --}}}  
        --{{{  Fourier transform in time direction
        VAL n IS 1.0 (REAL32) / (REAL32 ROUND l.t) :
        SEQ t = 0 FOR l.t
          SEQ k.x = 0 FOR l.s
            SEQ k.y = 0 FOR l.s
              SEQ k.z = 0 FOR l.s
                SEQ
                  s.1 [t] [k.x] [k.y] [k.z] := 0.0 (REAL32)
                  s.t [t] [k.x] [k.y] [k.z] := 0.0 (REAL32)
                  s.x [t] [k.x] [k.y] [k.z] := 0.0 (REAL32)
                  s.y [t] [k.x] [k.y] [k.z] := 0.0 (REAL32)
                  s.z [t] [k.x] [k.y] [k.z] := 0.0 (REAL32)
                  SEQ k = 0 FOR l.t
                    VAL e IS re.phase [t] [k] :
                    VAL f IS im.phase [t] [k] :
                    SEQ
                      s.1 [t] [k.x] [k.y] [k.z] := s.1 [t] [k.x] [k.y] [k.z] +
                                              (e * r.1 [k] [k.x] [k.y] [k.z])
                      s.t [t] [k.x] [k.y] [k.z] := s.t [t] [k.x] [k.y] [k.z] -
                                              (f * r.t [k] [k.x] [k.y] [k.z])
                      s.x [t] [k.x] [k.y] [k.z] := s.x [t] [k.x] [k.y] [k.z] +
                                              (e * r.x [k] [k.x] [k.y] [k.z])
                      s.y [t] [k.x] [k.y] [k.z] := s.y [t] [k.x] [k.y] [k.z] +
                                              (e * r.y [k] [k.x] [k.y] [k.z])
                      s.z [t] [k.x] [k.y] [k.z] := s.z [t] [k.x] [k.y] [k.z] +
                                              (e * r.z [k] [k.x] [k.y] [k.z])
                  s.1 [t] [k.x] [k.y] [k.z] := n * s.1 [t] [k.x] [k.y] [k.z]
                  s.t [t] [k.x] [k.y] [k.z] := n * s.t [t] [k.x] [k.y] [k.z]
                  s.x [t] [k.x] [k.y] [k.z] := n * s.x [t] [k.x] [k.y] [k.z]
                  s.y [t] [k.x] [k.y] [k.z] := n * s.y [t] [k.x] [k.y] [k.z]
                  s.z [t] [k.x] [k.y] [k.z] := n * s.z [t] [k.x] [k.y] [k.z]
        --}}}  
    :
    --}}}  
    --{{{  PROC pion.propagator    ([] REAL32, VAL [] REAL32 ^ 5)
    PROC pion.propagator ([] REAL32 pion, VAL [] [] [] [] REAL32 s.1, s.t, s.x, s.y, s.z)
    
      -- Zero momentum pion propagator.
    
      VAL n IS (REAL32 ROUND (3 * 4)) / (REAL32 ROUND (l.s * (l.s * l.s))) :
    
      SEQ t = 0 FOR l.t
        SEQ
          pion [t] := 0.0 (REAL32)
          SEQ k.x = 0 FOR l.s
            SEQ k.y = 0 FOR l.s
              SEQ k.z = 0 FOR l.s
                VAL s.1 IS s.1 [t] [k.x] [k.y] [k.z] :
                VAL s.t IS s.t [t] [k.x] [k.y] [k.z] :
                VAL s.x IS s.x [t] [k.x] [k.y] [k.z] :
                VAL s.y IS s.y [t] [k.x] [k.y] [k.z] :
                VAL s.z IS s.z [t] [k.x] [k.y] [k.z] :
                pion [t] := pion [t] + ((s.1 * s.1) + ((s.t * s.t) +
                        ((s.x * s.x) + ((s.y * s.y) +  (s.z * s.z)))))
          pion [t] := n * pion [t]
    :
    --}}}  
    
    REAL32 k.fermion, r.fermion :
    INT r.ini, r.fin :  -- Smearing radii at initial and final points.
    [l.t] [l.s] [l.s] [l.s] REAL32 s.1, s.t, s.x, s.y, s.z :
    [l.t] REAL32 pion :
    SEQ
      --{{{  get parameters
      SEQ
        term.write (screen, "Pion propagator in free fermion theory*c*n*n")
        term.write.INT (screen, "Lattice size = {} ** ", l.t)
        term.write.INT (screen, "{}****3*c*n*n", l.s)
        term.read.REAL32 (keyboard, screen, "Hopping parameter?   : ", k.fermion)
        term.read.REAL32 (keyboard, screen, "r-parameter?         : ", r.fermion)
        term.read.INT    (keyboard, screen, "Initial smear radius?: ", r.ini)
        term.read.INT    (keyboard, screen, "Final   smear radius?: ", r.fin)
        term.write (screen, "*c*n")
      --}}}  
      fermion.propagator (s.1, s.t, s.x, s.y, s.z,
                          k.fermion, r.fermion, r.ini, r.fin)
      pion.propagator    (pion, s.1, s.t, s.x, s.y, s.z)
      --{{{  print pion
      INT any :
      SEQ
        SEQ i = 0 FOR l.t
          SEQ
            term.write.INT    (screen, "{2}: ", i)
            term.write.REAL32 (screen, "{E12.6}*c*n", pion [i])
        keyboard ? any
      --}}}  
    
    --}}}  
    term.close (screen)
  scrn.save (to.save, screen, from.user.filer [0], to.user.filer [0],
             "dump", "dump")
