//
//      Program sarep.pls
//
//      Plots sequential assignments
//
//      (C) Copyright 1996 by Carlsberg Laboratory
//          All rights reserved
//


//
//      Specify the area to analyse:
//

DIM chainid$ OF 8,tmp$ of 2
chainid$:=""
resmin%:=1
resmax%:=99
aaperline%:=33
nsp%:=1
jsp_id%:=20

DIM sslist%(resmin%:resmax%),sp_ids%(nsp%)
DIM alist$(resmin%:resmax%) OF 1
DIM se. OF se_rec
DIM aa. OF aa_rec

READ nlist%
DIM list.(nlist%) OF sa_setup
FOR i%:=1 TO nlist% DO
  READ list.(i%).delta%
  READ list.(i%).color%
  READ list.(i%).density
  READ list.(i%).showintensity%
  READ list.(i%).h1_name$
  READ list.(i%).h2_name$
NEXT i%


PAGE

sp_ids%(1):=21

//
//      Fill in the sslist% table:
//

FOR ires%:=resmin% TO resmax% DO
  sslist%(ires%):=-1
NEXT ires%

ss_id%:=CAT_SS_FIRST
WHILE ss_id%<>-1 DO
  ires%:=CAT_SEQ_POS_OF_SS(ss_id%)
  IF ires%>=resmin% AND ires%<=resmax% THEN
    IF chainid$="" OR CAT_SEQ_CHAIN_OF_SS(ss_id%)=chainid$ THEN
      sslist%(ires%):=ss_id%
    ENDIF
  ENDIF
  ss_id%:=CAT_SS_NEXT(ss_id%)
ENDWHILE

sc_id%:=CAT_SC_FIRST
IF chainid$<>"" THEN
  WHILE sc_id%<>-1 DO
    CAT_SC_GET sc_id%,sc.
    IF sc.sc_chainid!=chainid$ THEN
      GOTO scfound
    ENDIF
    sc_id%:=CAT_SC_NEXT(sc_id%)
  ENDWHILE
scfound:
ENDIF

IF sc_id%=-1 THEN
  PRINT "Could not find chain."
  STOP
ENDIF

nres%:=resmax%-resmin%+1

se_id%:=CAT_SC_SE_FIRST(sc_id%)
WHILE se_id%<>-1 DO
  CAT_SE_GET se_id%,se.
  IF se.se_position%>=resmin% AND se.se_position%<=resmax% THEN
    aa_id%:=CAT_AA_OF_SE(se_id%)
    CAT_AA_GET aa_id%,aa.
    alist$(se.se_position%):=aa.aa_name1!
  ENDIF
  se_id%:=CAT_SC_SE_NEXT(se_id%,sc_id%)
ENDWHILE

//
//      Start the plotting stuff:
//
PDEVICE "ps"
PINIT "sarep.ps",dev_xmin,dev_xmax,dev_ymin,dev_ymax
//
//      y:              Start y position at top
//      dy:             step between each line
//      dx:             step between each amino acid
//      charsize:       Character height
//
charsize:=0.3
y:=dev_ymax-charsize*2
dy:=charsize*1.5
dx:=charsize
xbase:=dev_xmin+12.5*charsize


FOR ires1%:=resmin% TO resmax% STEP aaperline% DO
  ires2a%:=ires1%
  ires2b%:=ires1%+aaperline%-1
  IF ires2b%>resmax% THEN ires2b%:=resmax%
  // Plot numbers&names
  PTEXTSTYLE charsize,4
  FOR ires2%:=ires2a% TO ires2b% DO
    IF ires2% MOD 5=0 OR ires2%=resmin% OR ires2%=resmax% THEN
      PTEXT xbase+(ires2%-ires2a%)*dx,y,STR$(ires2%)
    ENDIF
  NEXT ires2%
  y:=y-dy
  PTEXTSTYLE charsize,1
  PTEXT dev_xmin,y,"Sequence"
  PTEXTSTYLE charsize,4
  FOR ires2%:=ires2a% TO ires2b% DO
    PTEXT xbase+(ires2%-ires2a%)*dx,y,alist$(ires2%)
  NEXT ires2%
  y:=y-dy*1.5
  
  FOR k%:=1 TO nlist% DO
    PTEXTSTYLE charsize,1
    IF list.(k%).delta%=-1 THEN
      PTEXT dev_xmin,y,"J"
      FOR ires2%:=ires2a% TO ires2b% DO
        i%:=ires2%-ires2a%
        jv:=jvalue(sslist%(ires2%),jsp_id%)
        IF jv<>-1 THEN
          // Full box is 15 Hz...
          x0:=xbase+i%*dx
          y1:=jv/15
          IF y1>1 THEN y1:=1
          y1:=y1*dy
          PRECTANGLESTYLE list.(k%).color%,list.(k%).density
          PRECTANGLE x0,y,dx*0.75,y1
        ENDIF
      NEXT ires2%
    ELSE
      PTEXT dev_xmin,y,"d"+list.(k%).h1_name$+list.(k%).h2_name$+"(i,i+"+STR$(list.(k%).delta%)+")"
      FOR ires2%:=ires2a% TO ires2b% DO
        i%:=ires2%-ires2a%
        ires3a%:=ires2%
        ires3b%:=ires2%+list.(k%).delta%
        IF ires3b%<=resmax% THEN
          j%:=isconnected%(nsp%,sp_ids%(),sslist%(ires3a%),list.(k%).h1_name$,sslist%(ires3b%),list.(k%).h2_name$)
          IF j%<>-1 THEN
            CASE j% OF
            WHEN 0,3
              bwidth:=charsize
            WHEN 1
              bwidth:=charsize*0.33
            WHEN 2
              bwidth:=charsize*0.67
            ENDCASE 
            IF list.(k%).delta%<2 THEN
              x0:=xbase+i%*dx
              x1:=x0+dx
              y0:=0
            ELSE
              x0:=xbase+i%*dx
              x1:=x0+(list.(k%).delta%+1)*dx
              bwidth:=bwidth/(list.(k%).delta%*2+1)
              y0:=bwidth*(i% MOD (list.(k%).delta%+1))*2
            ENDIF
            y1:=bwidth+y0
            PRECTANGLESTYLE list.(k%).color%,list.(k%).density
            PRECTANGLE x0,y0+y,x1-x0,y1-y0
          ENDIF //effect found
        ENDIF //within limits
      NEXT ires2%
    ENDIF
    y:=y-dy
  NEXT k%
  y:=y-dy*0.5
NEXT ires1%
print "Press return to finish: ",
input tmp$
ppage 
pfinish

STRUCTURE sa_setup
  DIM delta%,color%,density,showintensity%
  DIM h1_name$ OF 16,h2_name$ OF 16
ENDSTRUCTURE

DATA 4
DATA -1,1,1,TRUE,"",""
DATA 1,1,1,TRUE,"HN","HN"
DATA 1,1,1,TRUE,"HA","HN"
DATA 1,1,1,TRUE,"HB","HN"

FUNC isconnected%(nsp%,REF sp_ids%(),ss_id1%,aname1$,ss_id2%,aname2$) CLOSED
  //
  //    Returns max. intensity class iff at least one effect between
  //    ss_id1% and ss_id2% exists in the spectrum list
  //    sp_ids%
  //
  //    Return -1 if no connection
  //
  //    For each atom in ss_id1%, locate connecting cross peaks.
  //    For each cross peak belonging to a sp_ids% spectrum, see
  //    if ss_id2% is connected.
  //
  DIM cp. OF cp_rec
  DIM sname1$ OF LEN(aname1$)
  DIM sname2$ OF LEN(aname2$)
  
  max_intensity%:=-1
  
  IF ss_id1%=-1 OR ss_id2%=-1 THEN RETURN -1
  
  ss_at_id1%:=CAT_SS_AT_FIRST(ss_id1%)
  WHILE ss_at_id1%<>-1 DO
    sname1$:=SS_AT_NAME$(ss_at_id1%)
    IF sname1$=aname1$ THEN
      at_id1%:=CAT_AT_OF_SS_AT(ss_at_id1%)
      //
      //        Get spin systems connected to this atom,
      //        excluding ss_id1%
      //
      at_cp_id1%:=CAT_AT_CP_FIRST(at_id1%)
      WHILE at_cp_id1%<>-1 DO
        cp_id%:=CAT_CP_OF_AT_CP(at_cp_id1%)
        sp_id%:=CAT_SP_OF_CP(cp_id%)
        FOR i%:=1 TO nsp% DO
          IF sp_id%=sp_ids%(i%) THEN
            GOTO thissp
          ENDIF
        NEXT i%
        GOTO skipcp1
thissp:
        CAT_CP_GET cp_id%,cp.
        at_cp_id2%:=CAT_CP_AT_FIRST(cp_id%)
        WHILE at_cp_id2%<>-1 DO
          at_id2%:=CAT_AT_OF_AT_CP(at_cp_id2%)
          IF at_id1%<>at_id2% THEN
            ss_at_id2%:=CAT_AT_SS_FIRST(at_id2%)
            WHILE ss_at_id2%<>-1 DO
              sname2$:=SS_AT_NAME$(ss_at_id2%)
              IF sname2$=aname2$ THEN
                IF CAT_SS_OF_SS_AT(ss_at_id2%)=ss_id2% THEN
                  IF cp.cp_intensity_class%>max_intensity% THEN
                    max_intensity%:=cp.cp_intensity_class%
                    IF max_intensity%>=3 THEN RETURN max_intensity%
                  ENDIF
                ENDIF
              ENDIF
              ss_at_id2%:=CAT_AT_SS_NEXT(ss_at_id2%)
            ENDWHILE
          ENDIF
          at_cp_id2%:=CAT_CP_AT_NEXT(at_cp_id2%)
        ENDWHILE
skipcp1:
        at_cp_id1%:=CAT_AT_CP_NEXT(at_cp_id1%)
      ENDWHILE
    ENDIF
    ss_at_id1%:=CAT_SS_AT_NEXT(ss_at_id1%)
  ENDWHILE //for each atom on ss_id1%
  RETURN max_intensity%
ENDFUNC isconnected%

FUNC jvalue(ss_id%,sp_id%) CLOSED
  // returns j-value for the HN-HA coupling
  // for a given spectrum and direction
  
  // return -1 if not found
  
  DIM ss. OF ss_rec
  DIM cp. OF cp_rec
  DIM sp. OF sp_rec
  DIM at_cp. OF at_cp_rec
  DIM sname1$ OF 2
  DIM sname2$ OF 2
  
  nj%:=0
  jsum:=0
  
  // Scan through the atoms of ss_id%
  // If atom name is HN, look through the
  // the cross peaks to get the peaks
  // belonging to the sp_id% spectrum
  // that connects to an HA peak in the
  // same spin system. Get the J-value!
  
  IF ss_id%=-1 THEN RETURN -1
  
  ss_at_id1%:=CAT_SS_AT_FIRST(ss_id%)
  WHILE ss_at_id1%<>-1 DO
    sname1$:=SS_AT_NAME$(ss_at_id1%)
    IF sname1$="HN" THEN
      at_id1%:=CAT_AT_OF_SS_AT(ss_at_id1%)
      at_cp_id1%:=CAT_AT_CP_FIRST(at_id1%)
      WHILE at_cp_id1%<>-1 DO
        cp_id%:=CAT_CP_OF_AT_CP(at_cp_id1%)
        IF sp_id%=CAT_SP_OF_CP(cp_id%) THEN
          CAT_CP_GET cp_id%,cp.
          at_cp_id2%:=CAT_CP_AT_FIRST(cp_id%)
          WHILE at_cp_id2%<>-1 DO
            at_id2%:=CAT_AT_OF_AT_CP(at_cp_id2%)
            IF at_id1%<>at_id2% THEN
              ss_at_id2%:=CAT_AT_SS_FIRST(at_id2%)
              WHILE ss_at_id2%<>-1 DO
                sname2$:=SS_AT_NAME$(ss_at_id2%)
                IF sname2$="HA" THEN
                  IF CAT_SS_OF_SS_AT(ss_at_id2%)=ss_id% THEN
                    CAT_AT_CP_GET at_cp_id1%,at_cp.
                    CASE at_cp.at_cp_type% OF
                    WHEN 0
                      jv:=cp.cp_j_w1#
                    WHEN 1
                      jv:=cp.cp_j_w2#
                    WHEN 2
                      jv:=cp.cp_j_w3#
                    WHEN 3
                      jv:=cp.cp_j_w4#
                    ENDCASE 
                    nj%:=nj%+1
                    jsum:=jsum+jv
                  ENDIF
                ENDIF
                ss_at_id2%:=CAT_AT_SS_NEXT(ss_at_id2%)
              ENDWHILE
            ENDIF
            at_cp_id2%:=CAT_CP_AT_NEXT(at_cp_id2%)
          ENDWHILE
        ENDIF
        at_cp_id1%:=CAT_AT_CP_NEXT(at_cp_id1%)
      ENDWHILE
    ENDIF
    ss_at_id1%:=CAT_SS_AT_NEXT(ss_at_id1%)
  ENDWHILE //for each atom on ss_id1%
  IF nj%>0 THEN
    RETURN jsum/nj%
  ELSE
    RETURN -1
  ENDIF
  
ENDFUNC jvalue