program M_TtoDST
c
c	   Obtiene Delta S_T de M ISOTERMA
C       via relación de Maxwell
c DST(Ti+Ti+1)/2,Bj) = 1/(Ti+1-Ti)*Suma(M(Ti+1,Bj)-M(Ti,Bj))*(Bj-Bj-1)
      character*50 fileM
	real*4 M(100,100), T(100),B(100,100),Ymag(100),Campo(100)
	real*4 Me, DST(20),Bc(20),MM(100),BB(100)
	integer NB(100)
	character*14 titulo
	Terr=0.25 ! Tolerancia en T
	Berr=0.01 !Tolerancia en B
	open(10,file='M_HtoDST.inp')
	read(10,*) rmasa
	read(10,*) NBc, iBup !iBup = 1 campo aumenta, -1 disminuye
	read(10,*) (Bc(k),k=1,NBc)
	do i=1,100
	read(10,*,end=1) T(i)
	read(10,'(A50)') fileM
	write(*,*) 'Leyendo fichero: ',fileM
	write(*,*) 'Temperatura: ',T(i)
	 open(12,file=fileM)
	write(*,*)
	write(*,*) fileM
	write(*,*)'============='
	pause
	titulo=' '
	goto 4
	if(abs(rmasa).lt.0.0001)then
	do while(titulo.ne.'INFO,HARMONICS')
	read(12,'(a14)') titulo
	enddo
	read(12,'(a14,f14.7)') titulo,rmasa
	rmasa=rmasa*1.e-3
	endif
	write(*,*)' masa(g): ',rmasa
c	pause
	do while(titulo.ne.'Comment,Time S')
	read(12,'(a14)') titulo !encabezamiento
	enddo
4     continue
	 do j=1,100
3     continue
c  datos del SQUID
	 read(12,*,err=3,end=2) basura,basura,Campo(j),Temp,Ymag(j)
c	write(*,*) Temp, Campo(j)
	if(abs(Temp-T(i)).gt.Terr) goto 3

c
c
c      Ojo que se fuerza M(B=0) =0
c
	if(abs(Campo(j)).lt.Berr*10000.)then
	Ymag(j)=0.
	Campo(j)=0.
	endif
      Campo(j)= Campo(j)*1.e-4  ! cambio de Oe a T
c
c  ficheros en 2 columas con B en Tesla
c	 read(12,*,err=3,end=2) B(i,j),M(i,j)
	 enddo
2	 close(12)
       NBT = j-1
	write(*,*) 'No puntos totales: ', NBT
c
c     elige los de  Bup o Ndown segun el caso
c
	jj=1
      do j=2,NBT
	DB=Campo(j)-Campo(j-1)
	if(DB*IBup.gt.0.) then
	if(jj.eq.1)then
	B(i,1)=Campo(1)
	jj=jj+1
	endif
	B(i,jj)=Campo(j)
	M(i,jj)=Ymag(j)
	jj=jj+1
	else
	cycle
	endif
	enddo !j
	NB(i)=jj-1
	do jj=1,NB(i)
	write(*,*) T(i),B(i,jj),M(i,jj)
	enddo
	write(*,*) 'No de puntos para T = ',T(i),Nb(i)

C
C Ordena los puntos por orden creciente de campo
c     
	write(*,*) ' Puntos ordenados por el campo'
	if(abs(B(i,1)).gt.abs(B(i,NB(i)))) then
       do j=1,NB(i)
	 BB(j) = B(i,NB(i)-j+1)
	 MM(j) = M(i,NB(i)-j+1)
	 enddo
	 do j= 1,NB(i)
	 B(i,j) = BB(j)
	 M(i,j) = MM(j)
	write(*,*) B(i,j), M(i,j)
	 enddo
	endif
      enddo 
1	close(10)
      NT=i-1
	write(*,*) ' No temperaturas ',NT
	pause
	open(11,file='M_HtoDST.dat')
	write(11,100)' T/B    '
100   format(a8,4x,$)
	do k=1,Nbc
	write(11,101) Bc(k)
	enddo 
101   format(4x,F5.2,4x,$)
      write(11,99)
99    format()
	do i=1, NT-1
	 write(11,102) (T(i+1)+T(i))/2.
c	write(*,*)'Temperatura: ',(T(i+1)+T(i))/2. 
102    format(F8.3,$)
	 do k=1,NBc
	  Suma1=0.
	  do j= 1,NB(i)-1
	  if(B(i,j+1).le.Bc(k)+Berr) then
c	write(*,*) B(i,j),B(i,j+1), M(i,j),M(i,j+1)
	  Suma1=Suma1+0.5*(M(i,j+1)+M(i,j))*(B(i,j+1)-B(i,j))
	  else
	  Me=(M(i,j)*(Bc(k)-B(i,j))+M(i,j+1)*(B(i,j+1)-Bc(k)))
	>     /(B(i,j+1)-B(i,j))
	  Suma1=Suma1+0.5*(Me+M(i,j))*(Bc(k)-B(i,j))
	  exit
	  endif
	  enddo
	  Suma2=0.
	  do j= 1,NB(i+1)-1
	  if(B(i+1,j+1).le.Bc(k)+Berr)then
c	write(*,*) B(i+1,j),B(i+1,j+1), M(i+1,j),M(i+1,j+1)
	  Suma2=Suma2+0.5*(M(i+1,j+1)+M(i+1,j))*(B(i+1,j+1)-B(i+1,j))
	  else
	  Me=(M(i+1,j)*(Bc(k)-B(i+1,j))+M(i+1,j+1)*(B(i+1,j+1)-Bc(k)))
	>     /(B(i+1,j+1)-B(i+1,j))
	  Suma2=Suma2+0.5*(Me+M(i+1,j))*(Bc(k)-B(i+1,j))
	  exit
	  endif
	  enddo
c	write(*,*) ' Campo: ', Bc(k)
c	  write(*,*) 'Suma 1:', Suma1,' no campos: ',NB(i)
c	  write(*,*) 'Suma 2:', Suma2,' no de campos: ',NB(i+1)
c        pause
	  DST(k)= -(Suma2-Suma1)/(T(i+1)-T(i))/rmasa
        write(11,103) DST(k)
103     format(2x,F11.3,$)
	 enddo
       write(11,99)
	enddo
	close(11)
      end