program MTtoDST
c
c Obtiene delta S_T de datos de magnetzación a campo cte
c
c Se supone que están medidos a las msismas temperaturas 
c o muy próximas para cada campo.
c
      character*30 filein,title
	real*4 B(50),T(100,50),M(100,50),DMdT(100,50),TD(100,50)
	real*4 DST(100,50),TM(100),DM(50)
	real*4 aa(100), bb(100)
	integer NT(50) ! no de temperaturas para cada campo
	open(10,file='MTtoDST.inp')
	read(10,*) rmasa !gramos de muestra
	read(10,*) T1,DTC,T2 !Tini, incre, Tfin calcular
      TM(1)=T1
      do k=2,100
	TM(k)=TM(k-1)+DTC
	if(TM(k).gt.T2) exit
	enddo !k
	NTC=k-1
	do j=1,50
	read(10,*,end=1) B(j),fac !fac multiplica los datos
	read(10,'(A30)') filein
	open(12,file=filein)
	read(12,'(A30)') title
	write(*,*) ' Campo: ',B(j)
	!pause
	do i=1,50
2	read(12,*,err=2,end=3) basura,basura,T(i,j),basura,
     >      basura,basura,M(i,j)
      write(*,*) i, T(i,j),M(i,j)
	M(i,j)=M(i,j)*fac
      enddo !i
3	NT(j)=i-1
      write(*,*) NT(j),' puntos en el fichero'
      !pause
      close(12)
	enddo !j
1     close(10)
      NB=j-1
c
c  Ordena los puntos por temperatura
      do j=1,NB
      if(T(NT(j),j).lt.T(1,j)) then
	do i=1,NT(j)
	aa(i)=T(i,j)
	bb(i)=M(i,j)
	enddo !i
	do i=1,NT(j)
	T(i,j) = aa(NT(j)-i+1)
	M(i,j) = bb(NT(j)-i+1)
	enddo !i
	endif
	write(*,*)
      write(*,*) ' datos ordenados por temperatura'
	do i=1,NT(j)
	write(*,*) T(i,j),M(i,j)
	enddo !i
	!pause
	enddo !j
C
c      goto 50 ! salta esta parte
c     Obtiene (dM/dT)_B por dos puntos
      do j=1,NB
	do i=1, NT(j)-1
	dMdT(i,j)=(M(i+1,j)-M(i,j))/(T(i+1,j)-T(i,j))
	TD(i,j)=(T(i,j)+T(i+1,j))/2. !Asigna la Tmedia 
	enddo !i
	write(*,*)'Campo: ',  B(j)
	write(*,*) '    Tmedia     dM/dT_B'
	do i=1,NT(j)-1
	write(*,*) TD(i,j),DMdT(i,j)
	enddo !i
	!pause
	enddo !j
c
c
c Integra dM/dT_B para obtener DS_T
c	
50    continue
       T0=305.
 	write(*,*) 'Comprobamos Interpolación a ',T0,' K'
	write(*,*) 'T lista, B, DMDT inter'

      do k= 1,NTC
c
c     Interpola dM/dT entre los puntos más próximos
c
      do j=1,NB
c Busca los tres puntos más próximos a TM(k)
      x0=TM(K)
	do i=1,NT(j)-1
	if(T(i,j).lt.x0.and.x0.le.T(i+1,j))then	!if1
	write(*,*) x0,' i= ',i,' j = ',j
      x1 = T(i,j)
	y1=M(i,j)
	x2=T(i+1,j)
	y2=M(i+1,j)
	if(i.eq.1) then	 !if2
	x3=T(i+2,j)
	y3=M(i+2,j)
	else 
	 x3=T(i-1,j)
	 y3=M(i-1,j)
	 if(abs(x3-x0).gt.abs(T(i+2,j)-x0)) then   !if3
	  x3=T(i+2,j)
	  y3=M(i+2,j)
       endif	!3
	endif !2
      write(*,'(7(x,f10.4))') x0,x1,y1,x2,y2,x3,y3
c
c calcula la derivada en T(k) por tres puntos
c
      DM(j)=(y2-y1)*(x2+x3-2*x0)/(x2-x1)/(x3-x1)
	>     - (y3-y2)*(x1+x2-2*x0)/(x3-x2)/(x3-x1)
c	if(TM(k).gt.T0-6.and.TM(k).lt.T0+6)then	!if4
	write(*,*) TM(K),B(j),DM(j)
c	write(*,*) TD(i,j),dmdT(i,j),TD(i+1,j),dmdt(i+1,j)
c	endif !4
	exit
	endif !1
	enddo !i Tmedida
	if(j.eq.1)then
	DST(k,1)=DM(1)*B(1)/rmasa
	else
	DST(k,j)=DST(k,j-1)+0.5*(DM(j)+DM(j-1))*(B(j)-B(j-1))/rmasa
	endif
	enddo !j campo
	if(TM(k).gt.T0-0.5.and.TM(k).lt.T0+0.5)then
	write(*,*) ' Derivadas interpoladas a ',T0,' K'
	suma=0
	do jj=1,NB
	suma=suma+DM(jj)
	write(*,*) TM(k),B(jj),DM(jj)
	enddo
	suma=B(NB)/NB*(suma-(DM(1)+DM(NB))/2.)
	write(*,*)' Integral: ',suma/rmasa
	endif    
	enddo !k T calculada
c
c
	open(11,file='MTtoDST.dat')
	write(11,'(2x,a4,x,$)') 'T/B'
	do j=1, NB
	write(11,'(2x,f10.4,$)') B(j)
	enddo
	write(11,'()')
	do k=1,NTC
	write(11,'(f8.3,$)') TM(k)
	 do j=1,NB
	 write(11,'(2x,f10.5,$)') -DST(k,j)
	enddo !j
	write(11,'()')
 	enddo !k
	close(11)
	end