(defun mk-cube (d)
  (with-open-file (of "cube27x8.fhfc" :direction :output :if-exists :overwrite)
    (loop for xbig in '(0 -1 1) do
      (loop for ybig in '(0 -1 1) do
        (loop for zbig in '(0 -1 1) 
              with wfn
              do
          (setf wfn 1)
          (loop for xsmall from 0 to 1 do
            (loop for ysmall from 0 to 1 do
              (loop for zsmall from 0 to 1 do
                (format of "19F ~2D   ~4D ~4D ~4D     0 0~%" wfn 
                  (+ (* d 2 xbig) (* d xsmall))
                  (+ (* d 2 ybig) (* d ysmall))
                  (+ (* d 2 zbig) (* d zsmall)))
                (incf wfn)))))))))


(defun mk-lin (d n)
  (with-open-file (of (format nil "linear-3by~D.fhfc" n) 
                     :direction :output :if-exists :overwrite)
    (loop for zbig in '(1 0 2) do
      (when (zerop (rem zbig 2)) (format of "~%"))
      (loop for i from 0 to (1- n) do
        (format of "19F ~2D   0  0  ~4D     1 0~%" (1+ i)
          (* d (+ (* n zbig) i)))))))

; вращение вектора вокруг произвольной оси
; v - вектор
; a - угол
; x - единичный вектор, лежащий на оси вращения
(defun rotv (v a x)
  (let* ((l (aref x 0))
        (m (aref x 1))
        (n (aref x 2))
        (l2 (* l l))
        (m2 (* m m))
        (n2 (* n n))
        (i (aref v 0))
        (j (aref v 1))
        (k (aref v 2))
        (cos (cos a))
        (sin (sin a)))
  (make-array 3 :initial-contents
    (list (+ (* i (+ l2 (* cos (+ m2 n2)))) (* j (- (* l m) (* l m cos) (* n sin))) (* k (- (* l n) (* l n cos) (* m -1 sin))))
          (+ (* i (- (* l m) (* l m cos) (* n -1 sin))) (* j (+ m2 (* (+ l2 n2) cos))) (* k (- (* m n) (* m n cos) (* l sin))))
          (+ (* i (- (* l n) (* l n cos) (* m sin)))    (* j (- (* m n) (* m n cos) (* l -1 sin)))  (* k (+ n2 (* (+ l2 m2) cos))))))))

(defun vec (x y z)
  (make-array 3 :initial-contents (list x y z)))

(defun multv (v d)
  (map 'vector #'* v (vec d d d)))

(defun addv (a b)
  (map 'vector #'+ a b))

(defun difv (a b)
  (map 'vector #'- a b))

(defun scalp (a b)
  (reduce #'+ (map 'list #'* a b)))

(defun fixort (v ortv) ; поправить вектор v, чтобы он стал более ортогонален
  (let ((i             ; вектору ortv
      (position 0.000001 (map 'list #'abs ortv) :test #'<))
      (sp (scalp v ortv)))
    (cond 
      ((eql i 0) 
        (vec (- (aref v 0) (/ sp (aref ortv 0))) (aref v 1) (aref v 2)))
      ((eql i 1) 
        (vec (aref v 0) (- (aref v 1) (/ (scalp v ortv) (aref ortv 1))) (aref v 2)))
      ((eql i 2) 
        (vec (aref v 0) (aref v 1) (- (aref v 2) (/ (scalp v ortv) (aref ortv 2))))))))

(defun mk-cf (n step cshift rcf alpha beta tilt xxcs yycs zzcs)
; формирование структуры СnF2n вдоль оси z
; step - шаг вдоль оси z
; cshift - смещение углерода от оси
; rcf - расстяние C-F
; alpha  - угол F-C-F (около 109 град)
; beta - угол винтовой оси (чуть больше 180)
; tilt - поворот вектора F-F вокруг оси CM (M - середина вектора F-F)
  (format t "~%Please arrange atom blocks manually~%")
  (with-open-file (of (format nil "twisted-~D.fhfc" n) 
                     :direction :output :if-exists :overwrite)
    (let ((alpha-rad/2 (* alpha 0.01745329278 0.5))
          (beta-rad (* beta   0.01745329278))
          (tilt-rad (* tilt   0.01745329278))
          (f1s (make-array 3)) ; три блока векторов для ЦГУ
          (f2s (make-array 3))
          xx yy)               ; списки осей тензоров экр
      (loop for i from 0 to (1- (* n 3)) 
            with angle = 0
             and     z = 0
            do
              (let* 
                ((cmax (rotv (vec 1 0 0) angle (vec 0 0 1)))
                 (cm (multv cmax cshift))
                 (cf1nt (rotv (vec rcf 0 0) (+ angle alpha-rad/2) (vec 0 0 1)))
                 (cf2nt (rotv (vec rcf 0 0) (- angle alpha-rad/2) (vec 0 0 1)))
                 (cf1t (rotv cf1nt tilt-rad cmax))
                 (cf2t (rotv cf2nt tilt-rad cmax))
                 (c (addv cm (vec 0 0 z)))
                 (f1 (addv (addv cm cf1t) (vec 0 0 z)))
                 (f2 (addv (addv cm cf2t) (vec 0 0 z)))
                 (wfn (rem i n))
                 (blki (floor i n)))
                (push f1 (aref f1s blki))
                (push f2 (aref f2s blki))
                (when (zerop blki) ; тензоры есть только у основных атомов
                  (push cmax xx)   ; а правильно ли это?!
                  (push (fixort (difv cf1t cf2t) cmax) yy))
              (incf angle beta-rad)
              (incf z step)))
      (loop for ten from 1
            for xxc in xx
            for yyc in yy
            do 
            (format of "princ ~D ~D ~D ~D   ~D   ~D ~D ~D   ~D ~D~%"
              ten (aref xxc 0) (aref xxc 1) (aref xxc 2) xxcs
                  (aref yyc 0) (aref yyc 1) (aref yyc 2) yycs zzcs))
      (format of "~%")
      (loop for i in '(1 0 2) do
        (loop for wfn2 from 0
              for f1 in (aref f1s i)
              for f2 in (aref f2s i)
              with xxc = xx
              and  yyc = yy
            do
;                (format of "13C ~2D   ~4D  ~4D  ~4D     1 0~%" (+ (* 2 n) (* 2 wfn))
;                  (aref c 0) (aref c 1) (aref c 2))
                (format of "19F ~2D   ~4D  ~4D  ~4D     ~D 0~%" (+ 1 (* 2 wfn2))
                  (aref f1 0) (aref f1 1) (aref f1 2) (1+ wfn2))
                (format of "19F ~2D   ~4D  ~4D  ~4D     ~D 0~%" (+ 2 (* 2 wfn2))
                  (aref f2 0) (aref f2 1) (aref f2 2) (1+ wfn2)))))  ))