;;;calling

; $ gsc -cc-options "-O2 -march=native" ising 
; $ time gsi ising -e "(main 1000)"

;;;end calling

;;;compatibility

;bigloo compat
;(define random-integer random)

;bigloo has no random real generator!!! (define (random-real n) (random ?))

;guile compat
;guile does not support srfi-27 enough (define random-integer random) (define random-real random:uniform)

;;;end compat 

(define (SPIN+) +1)
(define (SPIN-) -1)

;;; coupling constant
(define J 1) ;ferromagnet
;(define J -1) ;antiferromagnet

;;; inverse temperature times coupling constant
(define (BETA*J) 0.3)
;(define (BETA*J) 0.4)
;(define (BETA*J) 0.5)

(define (display-spin spin)
  (cond 
   ((= spin (SPIN+)) (display "+"))
   ((= spin (SPIN-)) (display "-"))
   (#t display "?")))

;;; sample length in one dimension
(define L 40)

;;; sample dimension
(define D 2)

;;; total number of spins
;;; if this is a power of 2
;;; then modulo can be fast
(define N (expt L D))

;round N down to power of 2
;(set! N (expt 2 (inexact->exact (floor (/ (log N) (log 2))))))

(define GRID (make-vector N (SPIN+)))

(define (display-spatial-sample grid)
  (let loop ((site 0))
    (cond 
     ((< site N)
      (display-spin (vector-ref grid site))
      (if (= (- L 1) (modulo site L)) (newline))
      (loop (+ 1 site)))))
  (if (not (= (expt L D) N)) (newline))
  (newline))

(define (display-sample grid)
  (for-each display-spin (vector->list grid))
  (newline))

(define (dimension-loop site kons)
  (let loop ((d 0) (ret '()))
    (cond
     ((< d D)
      (loop (+ 1 d) (kons d site ret)))
     (#t ret))))
  
;;;helical boundary conditions
(define (helical-neighbors site)
  (dimension-loop site (lambda (d site ret)
			 (cons (modulo (+ (expt L d) site) N)
			       (cons (modulo (- site (expt L d)) N)
				     ret)))))

(define (helical-plus-neighbors site)
  (dimension-loop site (lambda (d site ret)
			 (cons (modulo (+ (expt L d) site) N)
			       ret))))

;;; Metropolis => equal probabilities for each transition (single-spin-flip)
(define (random-site)
  (random-integer N))

(define (flip! grid site)
  (vector-set! grid site (* -1 (vector-ref grid site))))

(define (energy-jump/J grid site)
  (* 2 (vector-ref grid site)
     (apply + (map (lambda (n) (vector-ref grid n)) (helical-neighbors site)))))

;;; Metropolis => always accept energy-lowering transitions
(define (accept-flip grid site)
  (if (> 0 (energy-jump/J grid site))
      #t
      (< (random-real) (exp (* (BETA*J) (energy-jump/J grid site))))))

(define (single-spin-flip! grid)
  (let ((site (random-site)))
    (if (accept-flip grid site)
	(flip! grid site))))

(define (magnetization grid)
  (apply + (vector->list grid)))

(define (energy/J grid)
  (define (-energy/J-plus site)
    (* (vector-ref grid site)
       (apply + (map (lambda (n) (vector-ref grid n))
		     (helical-plus-neighbors site)))))
  (let loop ((s 0) (E 0))
    (cond
     ((< s N)
      (loop (+ s 1) (- E (-energy/J-plus s))))
     (#t
      E))))

(define (energy grid)
  (* J (energy/J grid)))

(define (time-step)
  (let loop ((n 0))
    (cond
     ((< n N)
      (single-spin-flip! GRID)
      (loop (+ n 1))))))
  
(define (display-params . p)
  (let loop ((p p))
    (cond
     ((not (null? p))
      (display (car p))(display ": ")(display (cadr p))
      (cond
       ((not (null? (cddr p)))		; not last elt?
	(display ", ")
	(loop (cddr p)))))))
  (newline))

(define (display-config)
  (display-params "D" D "N" N "L" L "BETA*J" (BETA*J) "RANDOM-STATE" (random-source-state-ref default-random-source)))

(define TIMES #f)

(define (main args)
  (random-source-pseudo-randomize! default-random-source 0 0)
;  (write args)(newline)
  (set! TIMES args)
  (display-config)
;  (display-spatial-sample GRID)
  (let loop ((t 0))
    (cond
     ((< t TIMES)
      (time-step)
      (loop (+ t 1)))))
  (display-params "E" (energy GRID) "M" (magnetization GRID))
  (display-spatial-sample GRID)
)
