subroutine chromatic_function
c evaluate chromatic function as defined in MAD manual
c    W=a**2+b**2
c  b=(1/beta)(d beta/d delta), a=(d alpha/d delta)-alpha/beta (d beta/d delta)
c
      IMPLICIT REAL*8(A-H,O-Z), INTEGER (I-N)
	include 'input.inc'
	character*72 newname
	include 'elements.inc'
	include 'tempstor.inc'
	include 'analc.inc'
      COMMON/CONST/PI,TWOPI,CRDEG,CMAGEN,CLIGHT,EMASS,ERAD,ECHG
      COMMON/FUNC/BETAX,ALPHAX,ETAX,ETAPX,BETAY,ALPHAY,ETAY,ETAPY,
     <  STARTE,ENDE,DELTAE,DNU,MFPRNT,KOD,NLUM,NINT,NBUNCH
      COMMON/TWF/BETAOX,ALPHOX,ETAOX,ETAPOX,ANUX,
     >            BETAOY,ALPHOY,ETAOY,ETAPOY,ANUY,IE
      COMMON/LUM/ UO,TAUX,ALPHX,ALPHY,TAUY,
     <  ALPHE,TAUE,SIGE,SIGX,EPSX,
     <SIGY,EPSY,SIGXT,SIGYT,ENERGY,TAUREV
	logical first
	data first/.true./
	if(first)then
	  call strip_dirname(filename,istart,iiend,leftbrak,ightbrak)
	  type 2731,filename(istart:iiend)//'_chrom'
2731	  format(1x,' Open ',a,
     &	  ' to write chromatic function ')
	  newname=filename(istart:iiend)//'_chrom'
	  open(unit=94,file=filename(istart:iiend)//'_chrom',type='new')
	  first=.false.
5102  WRITE(94,5120)
5120  FORMAT(/,' ELEMENT #    BETAX    BETAY   ',
     <   '   ETAX     ETAPX     ETAY     ETAPY      dbetax   dalphx',
     <   '  wx      dbetay   dalphy      wy')
c      WRITE(94,880)(NAME(NEL,IN),IN=1,4),IE,BETAX,ALPHAX,BETAY,
c     <ALPHAY,ETAX,ETAPX,ETAY,etapy,beta1x,alph1x,wx,beta1y,alph1y,wy
	endif

	nel=norlst(ie)
      GAMMAX = (1.0D0+ALPHAX*ALPHAX)/BETAX
      GAMMAY = (1.0D0+ALPH1Y*ALPH1Y)/BETAY
      GAMM1X = 2.d0*ALPH1X*ALPHAX/BETAX-gammax/betax*beta1x
      GAMM1Y = 2.d0*ALPH1Y*ALPHAY/BETAY-gammay/betay*beta1y
      XM1112 = TEMP(1,1)*TEMP(1,2)
      XM2111 = TEMP(2,1)*TEMP(1,1)
      XM1221 = TEMP(1,2)*TEMP(2,1)
      XM1222 = TEMP(1,2)*TEMP(2,2)
      dm11=temp(1,12)+2.d0*temp(1,7)*etaox+temp(1,8)*etapox
     >      +temp(1,9)*etaoy+temp(1,10)*etapoy
      dm12=temp(1,17)+temp(1,8)*etaox+2.d0*temp(1,13)*etapox
     >      +temp(1,14)*etaoy+temp(1,15)*etapoy
      dm21=temp(2,12)+2.d0*temp(2,7)*etaox+temp(2,8)*etapox
     >      +temp(2,9)*etaoy+temp(2,10)*etapoy
      dm22=temp(2,17)+temp(2,8)*etaox+2.d0*temp(2,13)*etapox
     >      +temp(2,14)*etaoy+temp(2,15)*etapoy
      dXM1112 = dm11*TEMP(1,2)+temp(1,1)*dm12
      dXM2111 = dm21*TEMP(1,1)+temp(2,1)*dm11
      dXM1221 = dm12*TEMP(2,1)+temp(1,2)*dm21
      dXM1222 = dm12*TEMP(2,2)+temp(1,2)*dm22
      dBETAOX = 2.d0*TEMP(1,1)*dm11*BETAX-2.0D0*dXM1112*ALPHAX+
     <     2.d0*dm12*TEMP(1,2)*GAMMAX+
     <     TEMP(1,1)*TEMP(1,1)*BETA1X-2.0D0*XM1112*ALPH1X+
     <     TEMP(1,2)*TEMP(1,2)*GAMM1X
      dALPHOX = -dXM2111*BETAX+2.0D0*dXM1221*ALPHAX-
     <     dXM1222*GAMMAX-XM2111*BETA1X+(1.0D0+2.0D0*XM1221)*ALPH1X-
     <     XM1222*GAMM1X
C
      YM1112 = TEMP(3,3)*TEMP(3,4)
      YM2111 = TEMP(4,3)*TEMP(3,3)
      YM1221 = TEMP(3,4)*TEMP(4,3)
      YM1222 = TEMP(3,4)*TEMP(4,4)
      dm11=temp(3,9)*etaox+temp(3,14)*etapox
     <      +2.d0*temp(3,18)*etaoy+temp(3,19)*etapoy+temp(3,21)
      dm12=temp(3,10)*etaox+temp(3,15)*etapox
     <      +temp(3,19)*etaoy+2.d0*temp(3,22)*etapoy+temp(3,24)
      dm21=temp(4,9)*etaox+temp(4,14)*etapox
     <      +2.d0*temp(4,18)*etaoy+temp(4,19)*etapoy+temp(4,21)
      dm22=temp(4,10)*etaox+temp(4,15)*etapox
     <      +temp(4,19)*etaoy+temp(4,22)*etapoy+temp(4,24)
      dYM1112 = dm11*TEMP(3,4)+temp(3,3)*dm12
      dYM2111 = dm21*TEMP(3,3)+temp(4,3)*dm11
      dYM1221 = dm12*TEMP(4,3)+temp(3,4)*dm21
      dYM1222 = dm12*TEMP(4,4)+temp(3,4)*dm22
      dBETAOY = 2.0D0*temp(3,3)*dm11*betay-2.d0*dYM1112*ALPHAX+
     <     2.d0*dm12*TEMP(3,4)*GAMMAY+
     <     TEMP(3,3)*TEMP(3,3)*BETA1Y-2.0D0*YM1112*ALPH1Y+
     <     TEMP(3,4)*TEMP(3,4)*GAMM1Y
      dALPHOY = -dYM2111*BETAY+2.0D0*dYM1221*ALPHAY-
     <     dYM1222*GAMMAY-YM2111*BETA1Y+(1.0D0+2.0D0*YM1221)*ALPH1Y-
     <     YM1222*GAMM1Y

      ax=dalphox-alphox/betaox*dbetaox
      bx=dbetaox/betaox
      wx=dsqrt(ax**2+bx**2)
      ay=dalphoy-alphox/betaox*dbetaox
      by=dbetaoy/betaoy
      wy=dsqrt(ay**2+by**2)
      WRITE(94,880)(NAME(NEL,IN),IN=1,4),IE,BETAOX,BETAOY,
     <ETAOX,ETAPOX,ETAOY,etapoy,dbetaox,dalphox,wx,
     <dbetaoy,dalphoy,wy
880   FORMAT(X,4A1,I5,X,F10.3,F10.4,F10.3,3F10.4,6E10.2)
      return
c
c1     write(94,4),ie
4     format(1x,' cosx>1. at ',i4)
c    2 COSY=.5D0*(TEMP(3,3)+TEMP(4,4))
c      IF(DABS(COSY).GE.1.0D0)GOTO 3
c      SIN=DSQRT(1-COSY*COSY)
c      SIN=DSIGN(SIN,TEMP(3,4))
c      RMU0Y=DATAN2(SIN,COSY)
c      IF(RMU0Y.LT.0.0D0)RMU0Y=TWOPI+RMU0Y
c      RNU0Y=RMU0Y/TWOPI
c      RMU1Y=-.5D0*(A317+A419)/SIN
c      BETA1Y=(A319-BETAOY*COSY*RMU1Y)/SIN
c      ALPH1Y=(.5D0*(A317-A419)-RMU1Y*ALPHOY*COSY)/SIN
cc
c      ay=alph1y-alphoy/betaoy*beta1y
c      by=beta1y/betaoy
c      wy=ay**2+by**2
cc
c3     write(94,4),ie
c5     format(1x,' cosy>1. at ',i4)
      end