C      ********************************************************************
C      * THE PROGRAM COMPUTE THE HOMOGENEOUS ABSORPTION SPECTRUM AND THE  *   
C      * HOMOGENEOUS AND INHOMOGENEOUS EMISION SPECTRUM USING THE McCAMBER* 
C      * THEORY.                                                          *
C      * THE PROGRAM COMPUTE THE PARAMETERS OF A NUMBER OF GAUSSIAN       *    
C      * FUNCTION THAT ARE FITTING AN ABSORBTION INHOMOGENEOUS SPECTRUM   *
C      ********************************************************************
	
	
	character resp*1 ,stat*1,sta*1,title*10  
	common /bl/ y(2000),yc(1000),xdist(2000)
	common /al/dx,xx,xinit
	common /kl/kinit,kfinal
	common /xl/xlamel  
	common /nn/n1,n2,n3,n4,n5,n6,n7,n8 

	dimension x(100)
	dimension a(100),b(100),c(100)
	dimension sim(1000)                          
	real intr
	open(unit=20,file='spectin.dat')
	open(unit=19,file='specout.dat')
	open(unit=16,file='xaux.dat')               
	open(unit=17,file='auxx.dat') 
	open(unit=15,file='homog.dat') 
	open(unit=14,file='coefih.dat') 
	open(unit=13,file='coefom.dat') 
	open(unit=12,file='compar.dat')           
	open(unit=11,file='INIT.dat')
	

	write(6,*) ' THE TITLE'
	READ(5,2333) title
2333    format(a10) 
	write(6,2126)
2126    format(' ',/,'do you want same information about the program?'/
     1  5x,'(y/n)')
	 READ (*, '(A)') CCC
	   IF (ccc .EQ. 'n') THEN
	     go to 2128
	   ELSE
	      end if


	WRITE(6,2118)
	write(6,2119)
2118    FORMAT(1X,'******************************************************
     1  **********'//
     1  14X,'PROGRAM  ---SPECTRUM----'//)
2119    format(5x,'THIS PROGRAM COMPUTE THE HOMOGENEOUS ABSORPTION '/
     1  2X,'CROSS SECTION SPECTRUM AND THE HOMOGENEOUS AND '/ 
     1  2X,'INHOMOGENEOUS EMISSION CROSS SECTION SPECTRUM USING'/
     1  2X,'THE MCCAMBER THEORY'/
     1  2X,'IT USES AS INPUT DATA THE INTENSITY ABSORPTIOB SPECTRUM  '/
     1  2X,'AND DATA ABOUT THE WAVE GUIDE  '///)
	write(6,*) ' PRESS ENTER TO CONTINUE'

		read(5,2129)stat
2129    format(a4)

	WRITE(6,2120)
2120    FORMAT(4X,'THE INPUT DATA FILES ARE:'/
     1  2X,'1)-SPECTIN.DAT-IT CONTAINS THE INPUT SPECTRUM INTENSITIES'/
     1  2X,'PROFILE, EQUALY SPACED AS POSITION'//)
	WRITE(6,2121)
2121    FORMAT(2X,'2)-XAUXX-IT CONTAINS ON EACH LINE THE PARAMETERS'/
     1  2X,'OF ONE GAUSSIAN CURVE IN THE FOLLOWING ORDER '/
     1  2X,'-AMPLITUDE'/ 
     1  1X,'-POSITION'/
     1  2x,'-INVERS OF THE FWHM'///)
    



		write(6,*) ' PRESS ENTER TO CONTINUE'

			read(5,2130)sta
2130    format(a4)

		WRITE(6,2125)
2125    FORMAT(4X,'THE OUTPUT DATA FILES ARE:'//
     1  5X,'I-SPECTOUT.DAT'/
     1  2X,'FIRST COLUMN- THE WAVELENGTH [um]'/
     1  2X,'SECONT COLUMN-THE DERIVED ABSOPTION/EMITION CROSS SECTION'/ 
     1  2X,'PROFILE '/
     1  2X,'THIRD COLUMN -THE SIMULATED ABSOPTION/EMISION INHOMOGENEOUS'/
     1  2X,'CROSS SECTION PROFILE '/
     1  2X,'FOURTH COLUMN -THE SIMULATED ABSOPTION/EMISION HOMOGENEOUS'/
     1  2X,'CROSS SECTION PROFILE '//
     1  5X,'II- FILE XAUX.DAT '/
     1  2X,'THE INITIAL PARAMETERS OF THE GAUSSIAN FUNCTION AS'/  
     1  2X,'COMPUTED BY THE PROGRAMS'//
     1  5X,'III-HOMOG.DAT'/
     1  2X,'THE HOMOGENEOUS CROSS SECTION PROFILE  '//
     1  5X,'IV-COEFIH.DAT-THE VALUES OF THE BEST FIT COEFFICIENTS  '/
     1  2X,'FOR THE INHOMOGENEOUS CROSS SECTION PROFILE'//
     1  5X,' V-COEFOM.DAT-THE VALUES OF THE BEST FIT COEFFICIENTS  '/
     1  2X,'FOR THE HOMOGENEOUS CROSS SECTION PROFILE'///)

     


		write(6,*) ' PRESS ENTER TO CONTINUE'

			read(5,2130)sta

	 WRITE(6,*) ' THE PROGRAM WAS MADE BY DAN GROBNIC,'
	 WRITE(6,*) 'ASSISTENT PROFESOR AT THE POLITEHNIC' 
	 WRITE(6,*) 'UNIVERSITY BUCHAREST ,OCTOBER 1994'

2128    continue
	write(6,2131)
2131    FORMAT(4X,' IF YOU ARE READY TO START THE PROGRAM PRESS'/
     1  2X,'ENTER, ELSE CTRL C'//)
			read(5,2130)RESP








	write(6,10)
10      format(' the number of experimental points')
	 read(5,*) n  
	write(6,*) ' the wavelenght increment[um] and initial value [um]'
	read(5,*) dx,xinit
	write(6,*) 'the fluorescence lifetime [ms]'
	read(5,*) tau
	write(6,*) 'the refractive index of the guide'
	read(5,*) xind

	write(6,11)
11      format(' the inf limit of the spectrum')
	read(5,*)kinit 
	write(6,12)
12      format(' the sup limit of the spectrum')
	read(5,*)kfinal
	write(6,*) 'the number of the lines to be fitted'
	read(5,*) nmax

	write(6,*) '  wvlenght  increment lifetime  ref.index '
	write(6,*) '   [um]       [um]      [ms] ' 

	write(6,*) xinit,dx,tau,xind

	ymin=100000.
	ymax=0.
	do i=1,n
	xdist(i)=xinit+i*dx
	read(20,*)dumm, y(i)
c       y(i)=-y(i)
	if (y(i).lt.ymin) then
	ymin=y(i)
	else
	end if
	enddo
	do  i=1,n
	y(i)=-ymin+y(i)
	enddo

	do i=1,n
	if (abs(y(i)).gt.ymax) then
	ymax=y(i)
	xmax=xdist(i)
	else
	end if
	enddo
	do i=1,n
	y(i)=y(i)/ymax
	enddo
	ymax=1.


c       compute the effective width of the spectrum 
	intr=0.
	
	do i=kinit,kfinal-1
	intr=intr+(xdist(i)*y(i)+y(i+1)*xdist(i+1))*dx/2
	enddo
	write(*,*) 'the surf of the spect',intr
	xlamef=intr/ymax
	cff=sqrt(4*log(2))
	xlamel=xlamef/cff
	xlamel1=xlamel
	write(6,*) ' the effective waveline width',xlamef,xlamel
	write(6,*)'confirm the value of EFFEC FWHM or change it',xlamel

	read(5,*) xlamel
	coff=8*3.14*3*xind*xind*xlamel1*tau
	sigpk=10*xmax*xmax*xmax/coff
	do i=1,n
	y(i)=y(i)*sigpk
	enddo
	write(6,*) 'COMPUTE THE ABSORPTION(1) / EMISSION (0) '
	WRITE(6,*) 'CROSS SECTION PROFILE'
	READ(5,*) INDSEC

	IF(INDSEC.EQ.1) THEN
	GO TO 1911
	ELSE
	ENDIF
c       calculation of the emision cross section (Mc Camber)
	alfac=66.4/1.38

	do i=1,n
	sim(i)=y(i)
	ylamd=xinit+i*dx
	cfc=alfac*(1/ylamd-1/xmax)
C        write(6,*)'***', cfc,ylamd,xmax
	y(i)=0.91*y(i)*exp(-cfc)
C        write(6,*) y(i),sim(i)  
	write(12,*)ylamd, y(i),sim(i) 
	enddo

1911    CONTINUE

	n=3*nmax
	ddx=(kfinal-kinit)*dx/nmax
	kd=int(ddx/dx)

	write (6,*) 'the initial values of the variabiles'
	write(6,*)' a(i),b(i),c(i)'

	write(6,*) ' the progam  computes the initial values(1)' 
	write(6,*) 'or will  read the data from the aux.dat file (0)' 
	read(5,*) iopt2
	if(iopt2.ne.1) then 
	
	write (6,*) 'the initial values of the variabiles'
	write(6,*)' a(i),b(i),c(i)'

	do i=1,nmax

	read(17,*)a(i),b(i),c(i)            
	 
	write(6,*)a(i),b(i),c(i) 
	c(i)=1./c(i) 
	
	end do
	go to 100 
	else
	endif

	ffi=0
	do i=1,nmax
	  k=int(i*ddx/dx)   
	ffi=ffi+y(k)
	enddo

	sqpi=sqrt(3.14)
	xlamlef=100*xlamel/99
	cfmul=intr/xlamlef/sqpi/ffi
	WRITE(6,*) '****THE COEFICIENT####',cfmul

	do i=1,nmax
	b(i)=xinit+i*ddx
	k=int(i*ddx/dx)
	if(abs(y(k)/ymax).lt.0.0001) then
	a(i)=ymax/5.
	else
	a(i)=y(k)*2/3
c        a(i)=y(k)*cfmul
	endif
	c(i)=2./(ddx)
c        c(i)=1./xlamlef
	write(6,*) 'the ',i,' gaussian function'
	write(6,*) a(i),b(i),c(i) 
	 
	end do

	do i=1,nmax
c        write(6,*)'***', a(i),b(i),1./c(i) 
	write(16,*) a(i),b(i),1./c(i)
	enddo
	close(unit=16)
	write(6,*) 'replace same of this values y(1)/n(0) ?'
	read(5,*) iopt 

	if(iopt.eq.0)  go to 100

	do i=1,nmax
	write(6,*) 'if you want to change gaussian',i,'push 1'
	read(5,*) iopt1
	if(iopt1.eq.1) then
	write(6,*) 'a(i),b(i),c(i)'
	read(5,*) a(i),b(i),c(i) 
	else
	endif
	enddo

	do i=1,nmax
	write(6,*) 'the ',i,' gaussian function'
	write(6,*) a(i),b(i),c(i)
	end do
100     continue

c       THE INITIAL SPECTRUM APROXIMATION
	do i=1,nmax
	write(6,*) 'a,b,c',a(i),b(i),c(i)
	x(i)=a(i)
	x(nmax+i)=b(i)
	x(2*nmax+i)=c(i)
	end do

	write(6,*) 'THE INITIAL SPECTRUM APROXIMATION',n
	
	do i=kinit,kfinal         
	xx=xinit+i*dx
	  call spect(n,x,spc)
	
	write(11,*)xx, y(i),spc
	end do



	do i=1,nmax
	xl=1./c(i)
	if(xl.lt.xlamel) then
	c(i)=99./xlamel/100.
	write(6,*) 'c(i)',c(i)
	else
	endif
	enddo



	do i=1,nmax
	write(6,*) 'a,b,c',a(i),b(i),c(i)
	x(i)=a(i)
	x(nmax+i)=b(i)
	x(2*nmax+i)=c(i)
	end do

	do i=1,3*nmax
	write(6,*) 'x(i)',i,x(i)
	enddo

	write(6,*) 'give the order of variable to be kept const.'
	write(6,*) 'else push 0'
	write(6,*) 'n1'
	read(5,*) n1
	write(6,*) 'n2'
	read(5,*) n2
	write(6,*) 'n3'
	read(5,*) n3
	write(6,*) 'n3'
	read(5,*) n4
	write(6,*) 'n5'
	read(5,*) n5
	write(6,*) 'n6'
	read(5,*) n6
	n7=0
	n8=0





	write(6,*) '#####13 ####', x(13)
	read(5,*) isx
	call hookj(x,f,n)
	write(6,*) '#####13 ####', x(13)
	read(5,*) isx


	write(6,*) '****THE RESULTS***'

	do i=kinit,kfinal
	xx=xinit+i*dx
	  call spect(n,x,spc)
	  sim(i)=spc
c        write(19,*)xx, y(i),spc
	end do

	WRITE(14,2333) title 
	write(6,*) 'amplitude,position,SWFA'
	do i=1,nmax
	write(6,*) 'the gaussian',i
	a(i)=x(i)
	b(i)=x(nmax+i)
	c(i)=x(2*nmax+i)
	write(6,*) a(i),b(i),c(i)
	end do

	write(14,*) 'amplitude center  swfa'
	do i=1,nmax
c        write(6,*) '??????',a(i),b(i),1/c(i) 
	write(14,*) a(i),b(i),1/c(i)
	enddo

c       compute the omogenous cross section

C        cmax=0.
C        do i=1,nmax
C        if(cmax.lt.c(i))then
C        cmax=c(i)
C        else
C        endif
C        enddo
C        xlamm=1./cmax
C        if(xlamm.lt.xlamel) then
C        write(6,*) 'the gaus fitt. function <inhomogenous broadening'
C        go to 112
C        else 
C        endif 


	do i=1,nmax
	xl=1./c(i)
	if(xl.lt.xlamel) then

	write(6,*) 'the gaus fitt. function <inhomogenous broadening'
	go to 112

	
	write(6,*) 'c(i)',c(i)
	READ(5,*) INXX
	else
	endif
	enddo


	do i=1,nmax
	xlamm=1./c(i)
	coef=sqrt(xlamm*xlamm-xlamel*xlamel)
	a(i)=a(i)*xlamm/coef
	c(i)=1./coef
	
	enddo
	write(13,2333) title 
	write(13,*)  'the coeficients of omogenous broadening'
	do i=1,nmax
	write(13,*) 'a,b,c',a(i),b(i),1./c(i)
	x(i)=a(i)
	x(nmax+i)=b(i)
	x(2*nmax+i)=c(i)
	end do
	 write(15,2333) title 
	do i=kinit,kfinal
	xx=xinit+i*dx
	call spect(n,x,sp)
	
	write(15,*) xx,sp
	write(19,*) xx,y(i),sim(i),sp     
	end do

112     continue
	stop
	end

	subroutine fun(n,x,f)

	common /bl/ y(2000),yc(1000),xdist(2000)
	common /al/dx,xx,xinit            
	common /kl/kinit,kfinal
	common /xl/xlamel
	dimension x(100) 
c        dimension a(100),b(100),c(100)

	n3=n/3
c        write(6,*)'####',n3,n,xlamel
	do i=2*n3,n
	if(abs(x(i)).lt.0.0000001) then
	write(6,*)'x(i)<0',i,x(i)
	go to 343
	else
	endif
	xl=1./x(i)
c        write(6,*) i,xl,xlamel
	if(xl.lt.xlamel) then
	x(i)=99./xlamel/100.
	write(6,*) 'c(i)',x(i)
	else
	endif
	enddo
	
343    continue       
	do i=kinit,kfinal
	xx=xinit+i*dx
	
	call spect(n,x,sp)
c        write(6,*) 'n,x(1),sp' ,n,x(1),sp  
	yc(i)=sp
	
	end do

c        do i=kinit,kfinal,50 
c        write(6,*)'#####', y(i),yc(i)
c        end do

	st=0.
	sto=0.
	nto=0
	do k=kinit,kfinal
	tr=y(k)
	tp=yc(k)
c        tr=tr*tr
c        tp=tp*tp
	xx=tr-tp
	yy=(tr+tp)
	yy=(tr+tp)/2.
	xx=xx*xx
	st=st+xx*yy

c        yy=yy*yy
	nto=nto+1
	sto=sto+yy
	end do
	st=st/sto 
	st=sqrt(st)
	sto=sto/nto
	st=st/sto 
	f=st
c        fq=sqrt(f)
	

	nmax=n/3
	write(6,*)' THE ERROR AND THE VARIABLE ' ,f
c        do i=1,nmax
c        write(6,*) 'the gaussian',i
c        a(i)=x(i)
c        b(i)=x(nmax+i)
c        c(i)=x(2*nmax+i)
c        write(6,*) a(i),b(i),c(i)
c        end do

	
	return
	
	end


	
	subroutine hookj(x,f,n)

	dimension x(100),xb(100),x0(100),xx0(100),dx(100)
	dimension xbar(100),sxz(100),dxt(100)
	dimension xxy(100)
	common /nn/n1,n2,n3,n4,n5,n6,n7,n8
	write(6,*) '1#####13 ####', x(13)
      

	eps=0.05
	do k=1,n
	dxt(k)=x(k)/10.
	end do
	do k=1,n
	sxz(k)=x(k)
	x0(k)=sxz(k)
	x(k)=x0(k)
	xx0(k)=x0(k)
	dx(k)=dxt(k)
	xb(k)=x(k)
	xxy(k)=x(k)
	end do 
	call fun(n,x,f)
	ff=f
120     j=1
201     format(' al the ',i4,'iter the v of f is',f10.4,'and dx',f10.4)
130     do 150 i=1,n

	krt=(i-n1)*(i-n2)*(i-n3)*(i-n4)*(i-n5)*(i-n6)*(i-n7)*(i-n8)
	if(krt.eq.0) go to 150
	x(i)=x(i)+dx(i)
	call fun (n,x,fc)

	if(fc.lt.f) go to 140
	dxd=2*dx(i)
	ixt=1
	alf1=1  
133     continue
	if (ixt.gt.100) then
	x(i)=0.
	write(6,888)
888     format(' WARNING one of the variable was impose 0')
	go to 555
	else
	endif

	if(abs(alf1).lt.0.0000001) then
	write(6,*)'alf<0',alf1
	else
	endif

	x11=x(i)-dxd/alf1
	if (x11.lt.0.) then
	ixt=ixt+1
	alf1=alf1*2
	go to 133
	else
	x(i)=x11
	endif
	call fun(n,x,fc)

	if(fc.lt.f) go to 140
	x(i)=x(i)+dx(i)
	go to 150
140     f=fc
150     continue
      

	if(f.lt.ff) go to 190
	if(j.ne.1) go to 170
	if (abs(x(1)).lt.0.0001) go to 103 
	ddx=dx(1)/x(1)
	ddx=abs(ddx)

	if(ddx.le.eps) go to 210
	go to 105
103     continue        
	if (abs(dx(1)).le.abs(x(1)/50)) go to 210
105     continue
	do 160 k=1,n
	krt=(k-n1)*(k-n2)*(k-n3)*(k-n4)*(k-n5)*(k-n6)*(k-n7)*(k-n8)
	if(krt.eq.0) go to 160

	dx(k)=dx(k)/2
160     continue
     
     

170     do 180 k=1,n
	krt=(k-n1)*(k-n2)*(k-n3)*(k-n4)*(k-n5)*(k-n6)*(k-n7)*(k-n8)
	if(krt.eq.0) go to 180

	x(k)=xb(k)
180     continue
	go to 120
190     continue
	do 193 k=1,n
	krt=(k-n1)*(k-n2)*(k-n3)*(k-n4)*(k-n5)*(k-n6)*(k-n7)*(k-n8)
	if(krt.eq.0) go to 193

	xxy(k)=x(k)
	xbar(k)=xx0(k)
	xb(k)=x(k)
	bet=2.
556     x0(k)=bet*xb(k)-xbar(k)
	if (x0(k).lt.0.) then
	bet=bet*1.5
	go to 556
	else
	endif
	x(k)=x0(k)
	xx0(k)=xb(k)
	ff=f
193     continue
     
      

200     j=j+1
	it=it+1
	xit=it
	iq=1
	qq=abs(xit/10-int(xit/10))
c        qq=abs(xit/100-int(xit/100))
c        if (qq.lt. 0.01) then 
	if (qq.lt. 0.1) then
	do 295 k=1,n
295     write(6,296) k,x(k)
296     format(' x(',i2,')=',f10.5)
	write(6,217) f
	
	write(6,1010)
1010    format(' are you willing to go on ? Y/1/,N/0/')
	
	read(5,*) iq
	else 
	endif
	if(iq.eq.0) then 
	
	do i=1,n
	x(i)=xxy(i)
	end do
     
	go to 1015
	else
	endif
	if(it.gt.1000) go to 210
	go to 130               
210     continue
555     continue
1015    continue
	
	do 215 k=1,n
215     write(6,216) k,x(k)
216     format(' x(',i2,')=',f10.5)
	call fun(n,x,f)
	write(6,217) f
217     format(' the value of the hi2 function is',f10.5)
	write(6,333) it
333     format( ' number of iterations',i4)
	do i=1,n
	sxz(i)=x(i)
	end do
	return
	end


	subroutine spect(nn,x,spc)
	common /al/dx,xx,xinit            
	
	
	dimension x(100)
	dimension a(100),b(100),c(100)

	n=nn/3

	do i=1,n
	a(i)=x(i)
	b(i)=x(n+i)
	c(i)=x(2*n+i)
	end do

	spc=0.
	do i=1,n
	term=(xx-b(i))*c(i)
	term=-term*term
	term=a(i)*exp(term)
	spc=spc+term
c         write(6,*) 'xx,a,b,c trm',xx,a(i),b(i),c(i),term  
	end do
	
	return
	end
