PROGRAM RESTA
C       RESTA FUNCION SUAVE DE DATOS EXPERIMENTALES
      USE MSFLIB !Necesario para ejecutar programas externos
	LOGICAL(4) reswincom !resultado de ejecutar programas externos 
	CHARACTER*30 OLDFIL, FBASE
	CHARACTER*30 FICH
	COMMON /COEF/ NCOEF,A(30),TL,FS,NCOEFB,B(30),TLB,FSB,TCUT
	COMMON /THDEB/ TT(2000),THETA(2000),NTH
	COMMON /LANDAU/ BAT2,CAT3,TCR !Coeficientes de landau Y T (tri)crítica
      COMMON /SCHOTTKY/ NE,E(100),NG(100)
	COMMON /CBASE/ TBASE(20000),CPBASE(20000),NBASE
	DOUBLE PRECISION A,TL,FS,B,FSB,TLB
	DOUBLE PRECISION DEB,EIN,X
	DIMENSION TEIN(10)
	INTEGER NDEG(10),NRUN(2000)
	DIMENSION CP(20000),T(20000),DIF(20000),RES(20000)
	WRITE(6,*)' Resta funcion suave de datos experimentales de Cp'
	WRITE(6,*)'NUMERO DE COEFICIENTES(0 DE FICHERO)'
	WRITE(6,*)'-1 : FUNCION DE DEBYE'
	WRITE(6,*)'-2 : FUNCION DE EINSTEIN'
	WRITE(6,*)'-3 : NERNST-LINDEMANN; CP-CV = A*CP**2*T'
	WRITE(6,*)'-4 : LANDAU; C = 1/SQRT[4.*BB*BB + 12.*CC*(TC-T)]'
	WRITE(6,*)'-5 : Schottky'
      WRITE(6,*)'-6 : Resta de fichero sin suavizar base'
	READ(5,*) NCONTR
	NCOEF=NCONTR
	IF(NCONTR.GE.0) CALL LEECOEF
	IF(NCONTR.EQ.-1) GO TO 5
	IF(NCONTR.EQ.-2) GO TO 6
	IF(NCONTR.EQ.-3) GO TO 40
	IF(NCONTR.EQ.-4) GO TO 50
	IF(NCONTR.EQ.-5) GO TO 60
	IF(NCONTR.EQ.-6) GO TO 70
	GO TO 11
C
C        Datos de Theta de Debye (constante o variable con T,  seg'un tabla)
5       WRITE(6,*)' No de átomos por fórmula unidad Y THETA'
	WRITE(6,*)' Si THETA =< 0 se supone variable con T'
	READ(5,*) RDEB,THETA(1)
	NTH=1
	IF(THETA(1).LE.0.) THEN
	WRITE(6,*)' Fichero de datos de THETA(T) (dos columnas)'
	READ(5,100) FICH
	OPEN(10,file=FICH)
20      READ(10,*,ERR=20,END=21) TT(NTH),THETA(NTH)
	NTH=NTH+1       
	GO TO 20
21      CLOSE(10)
	NTH=NTH-1
	ELSE
	THD=THETA(1)
	ENDIF   
	GO TO 11
C
C       Datos para funci'on de Einstein
C
6       WRITE(6,*)' Número de modos INDEPENDIENTES'
	READ(5,*) NEIN
	WRITE(6,*) 'Temp. de Einstein y degeneración'
	READ(5,*) (TEIN(I),NDEG(I),I=1,NEIN)
	GO TO 11
C
C       Nernst-Lindemann
C
40      WRITE(6,*)'Constante A de Nernst- Lindemann'
	READ(5,*) ANERNST
	WRITE(6,*) ' Numero de coeficientes para calculo de Cp suave'
	WRITE(6,*) ' (0  para leer de fichero)'
	write(6,*) ' (-1 para tomar Cp sin suavizar)'
	READ(5,*) NCOEF
	IF(ncoef.ge.0) CALL LEECOEF
	GO TO 11
C
C       Landau
C
50      WRITE(6,*) 'Coeficientes de Landau: Tcr, B/At2, C/At3 '
	READ(5,*) TCR,BAT2,CAT3
	GOTO 11
C
C       Schottky
C
60    WRITE(6,*)'No de niveles'
 	READ(5,*) NE
	WRITE(6,*)' Energías en Kelvin y degeneraciones'
	READ(5,*) (E(I),NG(I), I=1,NE)
	WRITE(*,*)'No. de partículas por mol'
	READ(*,*) AP 
	GOTO 11
C
C  Resta datos de fichero sin suavizar
C  La y de la base para una x dada la obtiene interpolando entre los
C  datos más próximos
C
70    WRITE(6,*) 'Fichero datos de base para restar'
      READ(5,100) FBASE
	OPEN(13,FILE=FBASE)
	do nbase=1,20000
72	read(13,*,err=72,end=71) TBASE(NBASE),CPBASE(NBASE) 
	write(*,*)  TBASE(NBASE),CPBASE(NBASE)
	enddo
71	close(13)
      NBASE=NBASE-1
      write(*,*) 'Fin de lectura base: ', NBASE,' PUNTOS'
      GOTO 11
C
11      WRITE(6,*)'NOMBRE DE FICHERO DE DATOS'
	READ(5,100) OLDFIL      
100     FORMAT(A30)
	WRITE(6,*)'FORMATO ENTRADA: NRUN, T,Cp,... (0)'
	WRITE(6,*)'                  T, Cp,.....   (1)'
	READ(5,*) INP
	WRITE(6,*)'FORMATO SALIDA: NRUN, T,Cp,...(0)'
	WRITE(6,*)'                  T, Cp,.....   (1)'
	READ(5,*) IOU
	WRITE(6,*) 'Factor de escala: Cp/fs debe dar unidades de R'
	READ(5,*) FSCP
C
C       Estados correspondientes
C
	IF(NCONTR.GE.0) THEN
	WRITE(6,*)' Factor de escala en T, para multiplicar'
	READ(5,*) FMUL
	FS=FS*FMUL
	FSB=FSB*FMUL
	TL=TL*FMUL
	TLB=TLB*FMUL
	ENDIF
	OPEN(10,file=OLDFIL)
C	OPEN(11,file=NEWFIL)
	J=1
1       IF(INP.EQ.0) THEN
	READ(10,*,END=2,ERR=1) NRUN(J),T(J),CP(J)
	ELSE
	READ(10,*,END=2,ERR=1) T(J),CP(J)
	ENDIF
	CP(J)=CP(J)/FSCP
      IF(NCONTR.GE.0) RES(J)=F(T(J))
	IF(NCONTR.EQ.-3.AND.NCOEF.GE.0) RES(J)=ANERNST*F(T(J))**2*T(J)
	IF(NCONTR.EQ.-3.AND.NCOEF.LT.0) RES(J)=ANERNST*CP(J)**2*T(J)
	IF(NCONTR.EQ.-4) RES(J)=CLANDAU(T(J))
	IF(NCONTR.eq.-5) RES(J)=AP*CSH(T(J))
	IF(NCONTR.eq.-6) RES(J) = BASE(T(J))
	IF(NCONTR.EQ.-1) THEN
	   IF(NTH.GT.1) THD=THED(T(J))
	   IF(T(J).GT.0) THEN
	   X=DBLE(THD/T(J))
	   RES(J) = 3.*RDEB*SNGL(DEB(X))
	   ENDIF
	ENDIF
	IF(NCONTR.EQ.-2) THEN
	RES(J)=0.
	IF(T(J).GT.0.) THEN
	   DO 12 I=1,NEIN
	   X=DBLE(TEIN(I)/T(J))
12         RES(J)=RES(J)+NDEG(I)*SNGL(EIN(X))
	ENDIF
	ENDIF
	DIF(J)=CP(J)-RES(J)
	J=J+1
	GO TO 1
2       CLOSE(10)
	NJ=J-1

C
C       Escribe datos, resta y diferencia en RESTA.OUT
C
	OPEN(13,file='RESTA.OUT')
	IF(IOU.EQ.0) THEN 
	WRITE(13,*)' nrun      T(K)    	     Cp/R        Sust        Dif'
	ELSE
	WRITE(13,*)'       T(K)        Cp/R        Sust         Dif'
	ENDIF
	OPEN(17,FILE='RESTA.DAT')
	xmax=T(1)
	xmin=T(1)
	ymax=CP(1)
	ymin=CP(1)
      DO 14 J=1, NJ
	if(T(J).gt.xmax) xmax=T(J)
	if(T(J).lt.xmin) xmin=T(J)
	if(CP(J).gt.ymax) ymax=CP(J)
	if(CP(J).lt.ymin) ymin=CP(J)
	if(RES(J).gt.ymax) ymax=RES(J)
	if(RES(J).lt.ymin) ymin=RES(J)
	if(DIF(J).gt.ymax) ymax=DIF(J)
	if(DIF(J).lt.ymin) ymin=DIF(J)
	IF(IOU.EQ.0) THEN       
	WRITE(13,101) NRUN(J),T(J),CP(J),RES(J),DIF(J)
101	FORMAT(X,I4,4(2X,F10.3))
	ELSE
	WRITE(13,102) T(J),CP(J),RES(J),DIF(J)
102	FORMAT(4(2X,F10.3))
	ENDIF
	WRITE(17,*) T(J),DIF(J)
14      CONTINUE
      close(17)
	close(13)
c       genera RESTA.GNP y dibuja con gnuplot 

	OPEN(15,FILE='RESTA.GNP')
	write(15,*) 'set multiplot'
	write(15,*) 'set noautoscale'
	write(15,*) 'set xrange [',xmin,':',xmax,']'
	write(15,*) 'set yrange [',ymin,':',ymax,']'
	write(15,*) "plot '-'"
	do j=1,nj
	write(15,*) T(J),cp(J)
	enddo
	write(15,*) 'e'
	write(15,*) "plot '-' with lines"
	do j=1,nj
	write(15,*) T(J),RES(J)
	enddo
	write(15,*) 'e'
	write(15,*)"plot '-'"
	do j=1,nj
	write(15,*) T(J),dif(J)
	enddo
	write(15,*) 'e'
	CLOSE(15)
	reswincom=SYSTEMQQ('C:\calor\wgnuplot resta.gnp -')
	STOP
	END
C
C
	REAL FUNCTION F(T)
	COMMON /COEF/ NCOEF,A(30),TL,FS,NCOEFB,B(30),TLB,FSB,TCUT
	DOUBLE PRECISION A,TL,FS,B,TLB,FSB,DF,X
	IF(T.GT.TCUT.AND.TCUT.GT.0.) GO TO 5
	X=(T-TL)/FS
	IF(X.NE.0.D0) GO TO 2
	F=A(1)
	RETURN
2       DF=0.D0
	DO 1 I=1,NCOEF
1       DF=DF+A(I)*X**(I-1)
	F=SNGL(DF)
	RETURN
5       X=(T-TLB)/FSB
	IF(X.NE.0.D0) GO TO 4
	F=B(1)
	RETURN
4       DF=0.D0
	DO 3 I=1,NCOEFB
3       DF=DF+B(I)*X**(I-1)
	F=SNGL(DF)
	RETURN
	END
C
      DOUBLE PRECISION FUNCTION DEB(X)
      DOUBLE PRECISION DN,DD,X,ERR,XX
	DOUBLE PRECISION A(8)
C
C  LA EXPRESION RACIONAL DE THACHER (J. CHEM. PHYS. 32,638 ('60).
C  EL CALCULO ES VALIDO PARA VALORES DE 0<X<10, X=THETA/T. 
C       para x=12.5 da  el 1% de error
C
C       la segunda expresion es valida para x > 24, pero para x=12.5 da 
C       un 1% de error
C
	DEB=0.D0
	IF(X.LE.0.D0) RETURN
	IF(X.GT.15.D0) GO TO 1
	IF(X.GT.10.AND.X.LE.15) GO TO 3
      DN=3953.632D0-800.6087D0*X+85.07724D0*X**2-4.432582D0*X**3+
     >.0946173D0*X**4
      DD=3953.632D0+682.0012D0*X+143.155337D0*X**2+15.121491D0*X**3+
     >X**4
      DEB=4.D0*DN/DD-3.*X/(DEXP(X)-1.)
      RETURN
1     DEB=77.92727282/X**3
      IF (X.GT.30.) GO TO 2 
      ERR=3.*X/(DEXP(X)-1.)+DEXP(-X)*(12+36./X+72./X**2+72./X**3) 
      DEB=DEB-ERR 
2       RETURN
3       FS=3.72222D0
	TL=12.65000D0
	XX=(X-TL)/FS
      A(1)= 0.38325330510735D-01
      A(2)= -.33378319653339D-01
      A(3)= 0.18931752440859D-01
      A(4)= -.85545016616321D-02
      A(5)= 0.32348262820878D-02
      A(6)= -.10169656546186D-02
      A(7)= 0.23378056273646D-03
      A(8)= -.23474535567502D-04
	DEB=0.D0
	DO 4 I=1,8
4       DEB=DEB + A(I)* XX**(I-1)       
	RETURN
	END
C
C
      DOUBLE PRECISION FUNCTION EIN(X)
      DOUBLE PRECISION X
	EIN=0.D0
      IF (X.GT.30.D0) RETURN
      EIN=X**2*DEXP(X)/(DEXP(X)-1.)**2
	RETURN
	END
C
C
C
C
	SUBROUTINE INDEXX(N,ARRIN,INDX)
C
C       Esta subrutina ordena los indices
C       por valores crecientes de la variable ARRIN(I)
C       INDX = INDICES ORDENADOS: 
C       ARRIN(INDX(1)) < ARRIN(INDX(2)) < ARRIN(INDX(3) ...
C       Algoritmo: HEAPSORT. Numerical Recipes, p.233
C
	DIMENSION ARRIN(N),INDX(N)
	DO 11, J=1,N      !Inicializa INDX con numeros consecutivos
11      INDX(J)=J
	L=N/2+1
	IR=N               !Algoritmo HEAPSORT
10      CONTINUE
	IF(L.GT.1) THEN
	L=L-1
	INDXT=INDX(L)
	Q=ARRIN(INDXT)
	ELSE
	INDXT=INDX(IR)
	Q=ARRIN(INDXT)
	INDX(IR)=INDX(1)
	IR=IR-1
	IF(IR.EQ.1) THEN
	   INDX(1)=INDXT
	   RETURN
	   ENDIF
	ENDIF
	I=L
	J=L+L
20      IF(J.LE.IR) THEN
	   IF(J.LT.IR) THEN
		IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1))) J=J+1
	   ENDIF
	   IF(Q.LT.ARRIN(INDX(J))) THEN
		INDX(I)=INDX(J)
		I=J
		J=J+J
	   ELSE
		J=IR+1
	   ENDIF
	   GO TO 20
	ENDIF
	INDX(I)=INDXT
	GO TO 10
	END
C
C
	REAL FUNCTION THED(T)
C       Obtiene Theta de Debye a T, dada una tabla,
C        por interpolación (sin suavizados)
	COMMON /THDEB/ TT(2000),THETA(2000),NTH
	THED=THETA(NTH)
	IF(T.GE.TT(NTH)) RETURN
	THED=THETA(1)
	IF(T.LE.TT(1)) RETURN
	DO 1 I=2,NTH
	IF(T.LE.TT(I)) GO TO 2
1       CONTINUE
2       THED=THETA(I-1)+
     >       (T-TT(I-1))*(THETA(I)-THETA(I-1))/(TT(I)-TT(I-1))
	RETURN
	END

	SUBROUTINE LEECOEF
C       Lee/escribe coeficientes en fichero
	COMMON /COEF/ NCOEF,A(30),TL,FS,NCOEFB,B(30),TLB,FSB,TCUT
	DOUBLE PRECISION A,TL,FS,B,FSB,TLB
	CHARACTER*30 COEFIL
	IF(NCOEF.EQ.0) THEN
	WRITE(6,*)'Nombre de fichero de coeficientes'
	READ(5,100) COEFIL
100     FORMAT(A30)
	OPEN(12,file=COEFIL)
	READ(12,*) TCUT,NCOEF,FS,TL
	I=0
7       READ(12,*) A(I+1)
	I=I+1
	IF(I.LT.NCOEF) GO TO 7
	IF(TCUT.LE.0.) GO TO 8
	READ(12,*) NCOEFB,FSB,TLB
	I=0
70      READ(12,*) B(I+1)
	I=I+1
	IF(I.LT.NCOEFB) GO TO 70
8       CLOSE(12)
	ELSE
C
C       Datos de funci'on suave: coeficientes de polinomios ortogonales
C
4       WRITE(6,*)' Temperatura de corte (=0 si hay un solo tramo)'
	READ(5,*) TCUT
	DO 3 I=1,NCOEF
	IF(I.LT.10) WRITE(6,101) I-1
	IF(I.GE.10) WRITE(6,102) I-1        
101     FORMAT(' A(',I1,') =',$)
102     FORMAT(' A(',I2,') =',$)
	READ(5,*) A(I)
3       CONTINUE
	WRITE(6,*)'FACTOR DE ESCALA Y TRASLACION LATERAL'
	READ(5,*) FS,TL
	IF(TCUT.LE.0.) GO TO 31
	WRITE(6,*)' No. de coeficientes para T > Tcorte'
	READ(5,*) NCOEFB
	DO 30 I=1,NCOEFB
	IF(I.LE.10) WRITE(6, 201) I-1
	IF(I.GT.10) WRITE(6, 202) I-1        
201     FORMAT(' B(',I1,') =',$)
202     FORMAT(' B(',I2,') =',$)
	READ(5,*) B(I)
30      CONTINUE
	WRITE(6,*)'FACTOR DE ESCALA Y TRASLACION LATERAL'
	READ(5,*) FSB,TLB
C
C       guarda coeficientes en fichero
31      WRITE(6,*)'Nombre de fichero para guardar coeficientes'
	WRITE(6,*)' (RETURN = no guardar)'
	READ(5,100) COEFIL
	IF(COEFIL.EQ.' ') RETURN
	OPEN(15,file=COEFIL)
	WRITE(15,103) TCUT
103     FORMAT(1X,F7.2,T30,'Ta. de corte')
	WRITE(15,104) NCOEF
104     FORMAT(1X,I2,T30,'No. de coeficientes baja T')
	WRITE(15,105) FS
105     FORMAT(1X,D14.7,T30,'Factor de escala, baja T')
	WRITE(15,106) TL
106     FORMAT(1X,D14.7,T30,'Traslaci¢n lateral, baja T')
	DO 32 I=0,NCOEF-1
32      WRITE(15,107) A(I+1),I
107     FORMAT(1X,D14.7,T30,'A(',I2,')')
	WRITE(15,108) NCOEFB
108     FORMAT(1X,I2,T30,'No. de coeficientes alta T')
	WRITE(15,109) FSB
109     FORMAT(1X,D14.7,T30,'Factor de escala, alta T')
	WRITE(15,110) TLB
110     FORMAT(1X,D14.7,T30,'Traslaci¢n lateral, alta T')
	DO 33 I=0,NCOEF-1
33      WRITE(15,111) B(I+1),I
111     FORMAT(1X,D14.7,T30,'B(',I2,')')
	CLOSE(15)
	ENDIF
	RETURN
	END
	REAL FUNCTION CLANDAU(T)
C
C       Obtiene la capacidad calorífica an¢mala seg£n la teor¡a de Landau
C       para transiciones de 2o. orden o de 1o. d‚biles
C       La f¢rmula es 
C
C       (T/Cp)**2 = 4*BB*2 + 12.*CC*(TC-T) = 12.*CC*(TCR-T) (T<T0)
C
C       Con la energ¡a libre F = Fo + A(T-Tc)*q**2 + B*q**4 + C*q**6 + ...
C
C       y definiendo BB = B/A**2   y   CC = C/A**3
C
C       K.S.Aleksandrov y I.N. Flerov. SOV. PHYS. SOL.ST. 21(2) 195 (1979)
C
C       ATENCION !! La formula T0 - TC = 4*(TCR-T0) esta mal en Aleksandrov
C                Debe ser TCR-TC = 4*(TCR-T0) o T0-TC = 3*(TCR-T0)
C
C       Como consecuencia T0 = TCR - BB**2/12. (y no TCR - BB**2/15.)
	COMMON /LANDAU/ BB,CC,TCR !Coeficientes de landau Y T (tri)crítica
	IF(BB) 1,1,2
C       Transici¢n de 1er orden 
1       TC = TCR - BB*BB/CC/3.
	T0 = TCR - BB*BB/CC/12.
	GO TO 3
C
C       Transici¢n de 2o. orden
2       TC=TCR
	T0=TC
3       CLANDAU=0.
	IF(T.GT.T0) RETURN
	Y  = 4.*BB*BB + 12.*CC*(TC-T)
	CLANDAU = T/SQRT(Y)
	RETURN
	END
C
C
	REAL FUNCTION CSH(T)
      COMMON /SCHOTTKY/ NE,E(100),NG(100)
C
C       Calcula C de Schottky para niveles  E(I)
C       y degeneraciones NG(I)
C
	Z=0.
	U=0.
	U2=0.
	DO 1 I=1,NE
	Z=Z+NG(I)*EXP(-E(I)/T)
	U=U+E(I)*NG(I)*EXP(-E(I)/T)
1       U2=U2+E(I)*E(I)*NG(I)*EXP(-E(I)/T)
	CSH=(U2/Z-(U/Z)**2)/T/T
	RETURN
	END

	REAL FUNCTION BASE(T)
C
C Calcula el dato de base a una T dadA para restar
C 
	COMMON /CBASE/ TBASE(20000),CPBASE(20000),NBASE
	BASE=0.
	DO I=1,NBASE
	 if(T.eq.TBASE(I)) then 
	   BASE = CPBASE(I)
	   return
	 endif
	 IF((T.GT.TBASE(I).AND.T.LT.TBASE(I+1)).or.
	> (T.LT.TBASE(I).AND.T.GT.TBASE(I+1)))  then
	   t1=tbase(i)
	   t2=tbase(i+1)
	   BASE=(CPBASE(I)*(T2-T)+CPBASE(I+1)*(T-T1))/(T2-T1)
	 return 
	endif
	ENDDO
C   No hay datos de la base a esa T
	return
	end