//
//      helix1
//
//      Fits a cylinder to a set of coordinates.
//
PROC fithelix(st_id%,segment$,res1%,res2%,REF vector(),REF radius,REF length) CLOSED
  PROC vnorm(REF v())
    lng:=v(1)^2+v(2)^2+v(3)^2
    IF lng>0 THEN lng:=1/SQR(lng)
    v(1):=v(1)*lng
    v(2):=v(2)*lng
    v(3):=v(3)*lng
  ENDPROC vnorm
  PROC spherical(phi,theta,REF x,REF y,REF z)
    x:=SIN(theta)*COS(phi)
    y:=SIN(theta)*SIN(phi)
    z:=COS(theta)
  ENDPROC spherical
  FUNC opt4%(REF count%,cmax%,var%,REF xstart(),REF del0x(),REF xact(),REF xold(,),REF epsx(),yact,REF yold(),eps,mf) CLOSED
    PROC invert2(n%,REF a(,),eps) CLOSED
      DIM b(n%),c(n%),p%(n%),q%(n%)
      FOR k%:=1 TO n% DO
        pivot:=0
        FOR i%:=k% TO n% DO
          FOR j%:=k% TO n% DO
            IF ABS(a(i%,j%))>ABS(pivot) THEN
              pivot:=a(i%,j%)
              p%(k%):=i%
              q%(k%):=j%
            ENDIF
          NEXT j%
        NEXT i%
        IF ABS(pivot)<=eps THEN
          PRINT "pivot<=eps"
        ENDIF
        pk%:=p%(k%)
        IF pk%<>k% THEN
          FOR i%:=1 TO n% DO
            tmp:=a(k%,i%)
            a(k%,i%):=a(pk%,i%)
            a(pk%,i%):=tmp
          NEXT i%
        ENDIF
        qk%:=q%(k%)
        IF qk%<>k% THEN
          FOR i%:=1 TO n% DO
            tmp:=a(i%,k%)
            a(i%,k%):=a(i%,qk%)
            a(i%,qk%):=tmp
          NEXT i%
        ENDIF
        FOR j%:=1 TO n% DO
          IF j%=k% THEN
            b(j%):=1/pivot
            c(j%):=1
          ELSE
            b(j%):=-a(k%,j%)/pivot
            c(j%):=a(j%,k%)
          ENDIF
          a(k%,j%):=0
          a(j%,k%):=0
        NEXT j%
        FOR i%:=1 TO n% DO
          FOR j%:=1 TO n% DO
            a(i%,j%):=a(i%,j%)+c(i%)*b(j%)
          NEXT j%
        NEXT i%
      NEXT k%
      FOR k%:=n% TO 1 STEP -1 DO
        pk%:=p%(k%)
        IF pk%<>k% THEN
          FOR i%:=1 TO n% DO
            tmp:=a(i%,k%)
            a(i%,k%):=a(i%,pk%)
            a(i%,pk%):=tmp
          NEXT i%
        ENDIF
        qk%:=q%(k%)
        IF qk%<>k% THEN
          FOR i%:=1 TO n% DO
            tmp:=a(k%,i%)
            a(k%,i%):=a(qk%,i%)
            a(qk%,i%):=tmp
          NEXT i%
        ENDIF
      NEXT k%
    ENDPROC invert2
    PROC mult(n%,REF vec(),REF res())
      FOR i%:=1 TO n% DO
        r:=0
        FOR j%:=1 TO n% DO
          r:=r+matrix(i%,j%)*vec(j%)
        NEXT j%
        res(i%):=r
      NEXT i%
    ENDPROC mult
    PROC fit(dummy%)
      FOR i%:=1 TO fac% DO
        matrix(i%,1):=1
        k%:=2*var%+1
        FOR j%:=1 TO var% DO
          r:=xold(i%,j%)
          matrix(i%,j%+1):=r
          matrix(i%,1+var%+j%):=r*r
          IF j%<>var% THEN
            FOR m%:=j%+1 TO var% DO
              k%:=k%+1
              matrix(i%,k%):=r*xold(i%,m%)
            NEXT m%
          ENDIF
        NEXT j%
      NEXT i%
      invert2(fac%,matrix(,),eps)
      mult(fac%,yold(),coef())
    ENDPROC fit
    fac%:=(var%+1)*(var%+2) DIV 2
    DIM coef(fac%),vec(fac%),res(fac%)
    DIM matrix(fac%,fac%)
    corcount%:=count%-(count%-1) DIV fac%*fac%
    IF count%=0 THEN
      FOR i%:=1 TO var% DO
        xact(i%):=xstart(i%)
      NEXT i%
      GOTO out
    ENDIF
    yold(corcount%):=yact
    IF corcount%<fac% THEN
      FOR i%:=1 TO var% DO
        xact(i%):=xold(1,i%)
      NEXT i%
      i%:=corcount%
      r:=1
      IF corcount%<=2*var% THEN
        IF corcount%>var% THEN
          i%:=i%-var%
          r:=-1
        ENDIF
        xact(i%):=xact(i%)+r*del0x(i%)
      ELSE
        p%:=2*var%
        FOR i%:=1 TO var%-1 DO
          FOR j%:=i%+1 TO var% DO
            p%:=p%+1
            IF p%=corcount% THEN
              GOTO l2
            ENDIF
          NEXT j%
        NEXT i%
l2:
        xact(i%):=xact(i%)+del0x(i%)
        xact(j%):=xact(j%)+del0x(j%)
      ENDIF
      GOTO out
    ENDIF
    fit(1)
    FOR i%:=1 TO var% DO
      vec(i%):=-coef(i%+1)
      matrix(i%,i%):=2*coef(1+var%+i%)
    NEXT i%
    k%:=1+2*var%
    FOR i%:=1 TO var%-1 DO
      FOR j%:=i%+1 TO var% DO
        k%:=k%+1
        matrix(i%,j%):=coef(k%)
        matrix(j%,i%):=coef(k%)
      NEXT j%
    NEXT i%
    invert2(var%,matrix(,),eps)
    mult(var%,vec(),xact())
    FOR i%:=1 TO var% DO
      IF ABS(xact(i%)-xold(1,i%))>ABS(epsx(i%)) THEN
        GOTO out
      ENDIF
    NEXT i%
    RETURN 1
out:
    count%:=count%+1
    IF count%>cmax% THEN RETURN 2
    corcount%:=count%-(count%-1) DIV fac%*fac%
    IF corcount%=1 AND count%>1 THEN
      s:=0
      FOR i%:=1 TO var% DO
        r:=ABS((xact(i%)-xold(1,i%))/del0x(i%))
        IF r>s THEN
          n%:=i%
          s:=r
        ENDIF
      NEXT i%
      IF s>mf THEN
        FOR i%:=1 TO var% DO
          xact(i%):=xold(1,i%)+(xact(i%)-xold(1,i%))/s*mf
        NEXT i%
      ENDIF
    ENDIF
    FOR i%:=1 TO var% DO
      xold(corcount%,i%):=xact(i%)
    NEXT i%
    RETURN 0
  ENDFUNC opt4%
  FUNC fhelix(REF x())
    CASE fixdir% OF
    WHEN 1
      fx:=fixvalue
      fy:=x(1)
      fz:=x(2)
    WHEN 2
      fx:=x(1)
      fy:=fixvalue
      fz:=x(2)
    WHEN 3
      fx:=x(1)
      fy:=x(2)
      fz:=fixvalue
    ENDCASE 
    spherical(x(3),x(4),fdx,fdy,fdz)
    fr:=x(5)
    fsum:=0
    FOR fi%:=res1% TO res2% DO
      ft:=(c_ca(fi%,1)-fx)*fdx+(c_ca(fi%,2)-fy)*fdy+(c_ca(fi%,3)-fz)*fdz
      fsum:=fsum+(SQR((c_ca(fi%,1)-fx-ft*fdx)^2+(c_ca(fi%,2)-fy-ft*fdy)^2+(c_ca(fi%,3)-fz-ft*fdz)^2)-x(5))^2
      ft:=(c_c(fi%,1)-fx)*fdx+(c_c(fi%,2)-fy)*fdy+(c_c(fi%,3)-fz)*fdz
      fsum:=fsum+(SQR((c_c(fi%,1)-fx-ft*fdx)^2+(c_c(fi%,2)-fy-ft*fdy)^2+(c_c(fi%,3)-fz-ft*fdz)^2)-x(5))^2
      ft:=(c_n(fi%,1)-fx)*fdx+(c_n(fi%,2)-fy)*fdy+(c_n(fi%,3)-fz)*fdz
      fsum:=fsum+(SQR((c_n(fi%,1)-fx-ft*fdx)^2+(c_n(fi%,2)-fy-ft*fdy)^2+(c_n(fi%,3)-fz-ft*fdz)^2)-x(5))^2
      ft:=(c_hn(fi%,1)-fx)*fdx+(c_hn(fi%,2)-fy)*fdy+(c_hn(fi%,3)-fz)*fdz
      fsum:=fsum+(SQR((c_hn(fi%,1)-fx-ft*fdx)^2+(c_hn(fi%,2)-fy-ft*fdy)^2+(c_hn(fi%,3)-fz-ft*fdz)^2)-x(5))^2
    NEXT fi%
    RETURN fsum
  ENDFUNC fhelix
  DIM atc. OF atc_rec
  DIM c_ca(res1%:res2%,3)
  DIM c_c(res1%:res2%,3)
  DIM c_n(res1%:res2%,3)
  DIM c_hn(res1%:res2%,3)
  DIM v1(3),nv1(3),xact(5),xstart(5),del0x(5),epsx(5),xold(6*7 DIV 2,5),yold(6*7 DIV 2)
  FOR i%:=res1% TO res2% DO
    FOR j%:=1 TO 3 DO
      c_ca(i%,j%):=-9999
      c_c(i%,j%):=-9999
      c_n(i%,j%):=-9999
      c_hn(i%,j%):=-9999
    NEXT j%
  NEXT i%
  //read the coordinates:
  nres%:=res2%-res1%+1
  atc_id%:=CAT_ST_ATC_FIRST(st_id%)
  WHILE atc_id%<>-1 DO
    CAT_ATC_GET atc_id%,atc.
    IF segment$="" OR segment$=atc.atc_segment! THEN
      IF atc.atc_resno%>=res1% AND atc.atc_resno%<=res2% THEN
        CASE atc.atc_aname! OF
        WHEN "CA"
          c_ca(atc.atc_resno%,1):=atc.atc_x#
          c_ca(atc.atc_resno%,2):=atc.atc_y#
          c_ca(atc.atc_resno%,3):=atc.atc_z#
        WHEN "C"
          c_c(atc.atc_resno%,1):=atc.atc_x#
          c_c(atc.atc_resno%,2):=atc.atc_y#
          c_c(atc.atc_resno%,3):=atc.atc_z#
        WHEN "N"
          c_n(atc.atc_resno%,1):=atc.atc_x#
          c_n(atc.atc_resno%,2):=atc.atc_y#
          c_n(atc.atc_resno%,3):=atc.atc_z#
        WHEN "HN"
          c_hn(atc.atc_resno%,1):=atc.atc_x#
          c_hn(atc.atc_resno%,2):=atc.atc_y#
          c_hn(atc.atc_resno%,3):=atc.atc_z#
        OTHERWISE
        ENDCASE 
      ENDIF
    ENDIF
    atc_id%:=CAT_ST_ATC_NEXT#(atc_id%,st_id%)
  ENDWHILE
  //calculate initial guess:
  FOR i%:=1 TO 3 DO
    v1(i%):=0
    nv1(i%):=0
  NEXT i%
  FOR i%:=res1% TO res2% DO
    FOR j%:=1 TO 3 DO
      v1(j%):=v1(j%)+c_ca(i%,j%)+c_c(i%,j%)+c_n(i%,j%)+c_hn(i%,j%)
    NEXT j%
    IF i%>=res1%+4 THEN
      FOR j%:=1 TO 3 DO
        nv1(j%):=nv1(j%)+c_ca(i%,j%)-c_ca(i%-4,j%)
      NEXT j%
    ENDIF
  NEXT i%
  FOR j%:=1 TO 3 DO
    v1(j%):=v1(j%)/(4*nres%)
  NEXT j%
  vnorm(nv1())
  //print "center: ",v1(1);v1(2);v1(3)," direction: ",nv1(1);nv1(2);nv1(3)
  IF ABS(nv1(2))>ABS(nv1(1)) AND ABS(nv1(2))>ABS(nv1(3)) THEN
    fixdir%:=2
    xstart(1):=v1(1)
    fixvalue:=v1(2)
    xstart(2):=v1(3)
  ELIF ABS(nv1(3))>ABS(nv1(1)) AND ABS(nv1(3))>ABS(nv1(2)) THEN
    fixdir%:=3
    xstart(1):=v1(1)
    xstart(2):=v1(2)
    fixvalue:=v1(3)
  ELSE
    fixdir%:=1
    fixvalue:=v1(1)
    xstart(1):=v1(2)
    xstart(2):=v1(3)
  ENDIF
  theta0:=ACS(nv1(3))
  phi0:=ATN(nv1(2)/SIN(theta0),nv1(1)/SIN(theta0))
  //print "fixdir: ",fixdir%," theta0: ",theta0," phi0: ",phi0
  r0:=1.8345
  xstart(3):=phi0
  xstart(4):=theta0
  xstart(5):=r0
  // fit the helix:
  del0x(1):=0.01
  epsx(1):=0.0001
  del0x(2):=0.01
  epsx(2):=0.0001
  del0x(3):=0.01
  epsx(3):=0.0001
  del0x(4):=0.01
  epsx(4):=0.0001
  del0x(5):=0.01
  epsx(5):=0.0001
  
  count%:=0
h1:
  CASE opt4%(count%,5000,5,xstart(),del0x(),xact(),xold(,),epsx(),yact,yold(),1e-12,5000) OF
  WHEN 0
    //g2
    yact:=-fhelix(xact())
    //FOR i%:=1 TO 5 DO
    //PRINT xact(i%);
    //NEXT i%
    //PRINT yact
  WHEN 1
    //f1
    GOTO f1
  WHEN 2
    PRINT "Error!"
    GOTO f1
  ENDCASE 
  GOTO h1
  
f1:
  // calculate the extent of the helix:
  CASE fixdir% OF
  WHEN 1
    v1(1):=fixvalue
    v1(2):=xact(1)
    v1(3):=xact(2)
  WHEN 2
    v1(1):=xact(1)
    v1(2):=fixvalue
    v1(3):=xact(2)
  WHEN 3
    v1(1):=xact(1)
    v1(2):=xact(2)
    v1(3):=fixvalue
  ENDCASE 
  phi0:=xact(3)
  theta0:=xact(4)
  radius:=xact(5)
  spherical(phi0,theta0,vector(1),vector(2),vector(3))
  tlow:=1e+30
  thigh:=-tlow
  FOR i%:=res1% TO res2% DO
    t:=(c_ca(i%,1)-v1(1))*vector(1)+(c_ca(i%,2)-v1(2))*vector(2)+(c_ca(i%,3)-v1(3))*vector(3)
    IF t<tlow THEN tlow:=t
    IF t>thigh THEN thigh:=t
    t:=(c_c(i%,1)-v1(1))*vector(1)+(c_c(i%,2)-v1(2))*vector(2)+(c_c(i%,3)-v1(3))*vector(3)
    IF t<tlow THEN tlow:=t
    IF t>thigh THEN thigh:=t
    t:=(c_n(i%,1)-v1(1))*vector(1)+(c_n(i%,2)-v1(2))*vector(2)+(c_n(i%,3)-v1(3))*vector(3)
    IF t<tlow THEN tlow:=t
    IF t>thigh THEN thigh:=t
    t:=(c_hn(i%,1)-v1(1))*vector(1)+(c_hn(i%,2)-v1(2))*vector(2)+(c_hn(i%,3)-v1(3))*vector(3)
    IF t<tlow THEN tlow:=t
    IF t>thigh THEN thigh:=t
  NEXT i%
  length:=thigh-tlow
ENDPROC fithelix

DIM vector(3),vectors(4,3)
DIGITS 6
PAGE
OPEN FILE 1,"/indy2/mk/bk/helix/pracbp.dat",WRITE
PRINT FILE 1: "substr"+CHR$(9)+"v12"+CHR$(9)+"v13"+CHR$(9)+"v14"+CHR$(9)+"v23"+CHR$(9)+"v24"+CHR$(9)+"v34"

stc_id%:=1
st_id%:=CAT_STC_ST_FIRST(stc_id%)
st_no%:=1
WHILE st_id%<>-1 DO
  fithelix(st_id%,"",4,15,vector(),radius,length)
  FOR i%:=1 TO 3 DO
    vectors(1,i%):=vector(i%)
  NEXT i%
  PRINT st_id%;"Helix 1";vector(1);vector(2);vector(3);radius;length
  fithelix(st_id%,"",20,35,vector(),radius,length)
  FOR i%:=1 TO 3 DO
    vectors(2,i%):=vector(i%)
  NEXT i%
  PRINT st_id%;"Helix 2";vector(1);vector(2);vector(3);radius;length
  fithelix(st_id%,"",51,61,vector(),radius,length)
  FOR i%:=1 TO 3 DO
    vectors(3,i%):=vector(i%)
  NEXT i%
  PRINT st_id%;"Helix 3";vector(1);vector(2);vector(3);radius;length
  fithelix(st_id%,"",65,84,vector(),radius,length)
  FOR i%:=1 TO 3 DO
    vectors(4,i%):=vector(i%)
  NEXT i%
  PRINT st_id%;"Helix 4";vector(1);vector(2);vector(3);radius;length
  PRINT FILE 1: st_no%,CHR$(9),
  FOR i%:=1 TO 3 DO
    FOR j%:=i%+1 TO 4 DO
      PRINT FILE 1: ACS(vectors(i%,1)*vectors(j%,1)+vectors(i%,2)*vectors(j%,2)+vectors(i%,3)*vectors(j%,3))/PI*180,CHR$(9),
    NEXT j%
  NEXT i%
  PRINT FILE 1: 
  st_id%:=CAT_STC_ST_NEXT#(st_id%,stc_id%)
  st_no%:=st_no%+1
ENDWHILE
CLOSE FILE 1