* Copyright c 1998-2002 The Board of Trustees of the University of Illinois
* 		  All rights reserved.
* Developed by:	Large Scale Systems Research Laboratory
*               Professor Richard Braatz, Director*               Department of Chemical Engineering*		University of Illinois
*		http://brahms.scs.uiuc.edu
* * Permission hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to
* deal with the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the 
* Software is furnished to do so, subject to the following conditions:
* 		1. Redistributions of source code must retain the above copyright
*		   notice, this list of conditions and the following disclaimers.
*		2. Redistributions in binary form must reproduce the above 
*		   copyright notice, this list of conditions and the following 
*		   disclaimers in the documentation and/or other materials 
*		   provided with the distribution.
*		3. Neither the names of Large Scale Research Systems Laboratory,
*		   University of Illinois, nor the names of its contributors may
*		   be used to endorse or promote products derived from this 
*		   Software without specific prior written permission.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
* THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
* OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
*
*      param2.f
*
*	  This program reads measured concentration and moments data

*       from files and estimates crystal growth and nucleation 

*	  kinetic parameters using the weighted least-squares method.

*	  The data for this program are appropriate for the 

*	  batch crystallization of crystals with one characteristic

*       dimension.
*      
*       There are two input data files.  The first file, named 

*       concen_data11, should contain the concentration data.  In

*	  case you have second concentration data set, its name

*	  should be "concen_data12", and so on.	 The second file,

*	  named "mu_data11" should contain the measured moments data, 

*       in the format as follows:

*

*       moment00, moment10, moment01, moment11, moment20
*	  moment01, moment21, moment12, moment22, 
*	  moment30, moment31

*

*	  If you have a second set of moment data, then it should

*	  be named "mu_data12", and so on.

*

*	  If you are using this code for the ChE 391 course for a

*	  crystal with one or two characteristic dimensions, you 

*	  should only change CSTEP, MSTEP, Nsets, Cinterval,

*	  Minterval, x, weights, and Temp:

*

*	  CSTEP = total points in concentration measurements

*	  MSTEP = total points in moments measurements

*	  Nsets = numbers of data sets

*	  Cinterval = time interval between two concentration 

*				  measurements

*	  Minterval = time interval between two moment measurements

*	  x  =  initial guess of parameters

*	  weights = there are 6 weights to quantify the accuracy of

*			    the measurements.  Each weight is set to the

*			    inverse of the measurement error variance.

*				Normally it is very difficult to obtain an

*				accurate estimate of the zeroth order moment

*				and the fourth order moments, so the

*				corresponding weights for these moments are

*				normally set to zero.

*	  Temp = temperature profile (T(t))

*		  		

*	  Parameter changes	should be done both in the main 

*	  program and in the subroutine FCN.  Changes to Temp

*	  should be made in the subroutine Temp.

*

*	  If you are using this parameter estimation code for 

*	  another process, then you would need to change the model

*	  equations as well as the above variables.

*

*	  Output of the program is the best-fit model parameters,

*	  which is the vector called "x" in the output screen.  The
*       other parameters on the output screen give information on

*       the convergence of the parameter estimation algorithm,

*	  which can be ignored for most purposes.  The code takes

*	  up to 10 minutes to run, depending on the speed of the

*       computer and the quality of the initial parameter guesses. 

* 
*	The main program initializes the variables used by the
*	FFSQP subroutine, which minimizes the least-squared errors.  

*	Explanation of the variables is given in the FFSQP subroutine.

*	After initialization the program calls the FFSQP subroutine 

*	to solve the nonlinear constrained optimization problem.

*     Our experience is that the program converges provided that

*     you enter reasonable initial guesses.
*

*       Date:    February 12, 1998

*       Authors: Serena H. Chung and Richard D. Braatz

*                Department of Chemical Engineering

*                University of Illinois at Urbana-Champaign

*	  Modified:March 4, 2000

*	  By:      David L. Ma and Richard D. Braatz 

*
	use MSIMSL

      INTEGER CSTEP, MSTEP, Nsets, Cinterval, Minterval
      PARAMETER (CSTEP=160,Nsets=1, MSTEP=6)
      PARAMETER (Cinterval=1,Minterval=30)
      INTEGER nparam, npos1, npos2,I1
      PARAMETER (nparam = 6, npos1=12, npos2=8)
      CHARACTER*60 file_concen, file_mu, bob

      INTEGER  nf, nineqn, nineq, neqn, neq, iwsize, nwsize
      INTEGER mode, iprint,miter
      INTEGER inform, I
      PARAMETER(nf=1, nineq=0, neq=0)
      PARAMETER(mode=100,miter=10000)
      PARAMETER(iwsize=6*nparam+8*(nineq+neq)+7*(nf)+30) 
      PARAMETER(nwsize=4*nparam**2+5*(nineq+neq)*nparam+
     &		3*(nf)*nparam+
     &		26*(nparam+nf)+45*(nineq+neq)+100) 
      INTEGER iw(iwsize)
      REAL*8 bigbnd, eps, epseqn, udelta
      REAL*8 x(nparam), bl(nparam), bu(nparam)
      REAL*8 f(nf),g(nineq+neq+1),w(nwsize)
      INTEGER NWRITE, flag

      REAL*8 ConcData(Nsets,CSTEP), Mu00Data(Nsets,MSTEP)
      REAL*8 Mu10Data(Nsets,MSTEP), Mu01Data(Nsets,MSTEP)
      REAL*8 Mu11Data(Nsets,MSTEP), Mu20Data(Nsets,MSTEP)
      REAL*8 Mu02Data(Nsets,MSTEP), Mu21Data(Nsets,MSTEP)
      REAL*8 Mu12Data(Nsets,MSTEP), Mu22Data(Nsets,MSTEP)
	REAL*8 Mu30Data(Nsets,MSTEP), Mu31Data(Nsets,MSTEP)
	REAL*8 weight_conc, weight_mu00, weight_mu10, weight_mu01,
     &			 weight_mu11, weight_mu20, weight_mu02,
     &			 weight_mu21, weight_mu12, weight_mu22,
     &		     weight_mu30, weight_mu31 
 

      EXTERNAL FCN,cntr,grobfd,grcnfd

      COMMON /CONCEN/ConcData
      COMMON /MU00/Mu00Data
      COMMON /MU10/Mu10Data
      COMMON /MU01/Mu01Data
      COMMON /MU11/Mu11Data 
      COMMON /MU20/Mu20Data
	COMMON /MU02/Mu02Data
      COMMON /MU21/Mu21Data
      COMMON /MU12/Mu12Data 
      COMMON /MU22/Mu22Data
	COMMON /MU30/Mu30Data 
      COMMON /MU31/Mu31Data

      COMMON /DATA1/flag
      COMMON /DATA_OUT/NWRITE 

      COMMON /WEIGHT/weight_conc,
     &			 weight_mu00, weight_mu10, weight_mu01,
     &			 weight_mu11, weight_mu20, weight_mu02,
     &			 weight_mu21, weight_mu12, weight_mu22,
     &		     weight_mu30, weight_mu31 
 


*  Initial guess:

      DATA x/	1.15D0, 7.5D0, 1.65D0, 8.8D0,
     &		1.8D0, 18.174D0/


	print*, x(1),x(2),x(3),x(4),x(5),x(6)
	
* file name

      file_concen='concen_data'
      file_mu='mu_data'

*Input the weights for least-squares parameter estimation


      weight_conc=1.0D0/0.49D0/0.49D0
      weight_mu00=0.0D0
      weight_mu10=1.0D0/59008.0D0/59008.0D0
      weight_mu01=1.0D0/199721.0D0/199721.0D0
	weight_mu11=1.0D0/44892500.0D0/44892500.0D0
      weight_mu20=1.0D0/15474300.0D0/15474300.0D0
	weight_mu02=1.0D0/137020000.0D0/137020000.0D0
	weight_mu21=0.0D0
	weight_mu12=0.0D0
	weight_mu22=0.0D0
	weight_mu30=0.0D0
	weight_mu31=0.0D0
*FFSQP parameters


      bigbnd=1.0D12
      eps=1.0D-8
      epseqn=0.0D0
      udelta=0.0D0
      iprint=2
      nineqn=0
      neqn=0

* Lower bounds:


  	bl(1)=1.0D0
      bl(2)=5.0D0
      bl(3)=1.0D0
      bl(4)=5.0D0
	bl(5)=1.0D0
	bl(6)=3.0D0

* Upper bounds:


  	bu(1)=2.0D0
      bu(2)=15.0D0
      bu(3)=2.0D0
      bu(4)=15.0D0
	bu(5)=3.0D0
	bu(6)=20.0D0

	DO I =1+10 , Nsets+10

	    I1=I-10
* Read in data from file
        DO  J = 1,CSTEP
		   WRITE(file_concen(npos1:npos1+1), '(I2)') I
		   file_concen = file_concen(1:npos1+1) // '.txt'
             OPEN(UNIT=20, FILE=file_concen, FORM='FORMATTED',
     &             ACCESS='SEQUENTIAL', STATUS='OLD')
             READ(20,*) ConcData(I1,J)
	  ENDDO
	  CLOSE(UNIT=20)
	 
	  
	  DO J=1,MSTEP
		   WRITE(file_mu(npos2:npos2+1), '(I2)') I
		   file_mu = file_mu(1:npos2+1) // '.txt'
             OPEN(UNIT=20, FILE=file_mu, FORM='FORMATTED',
     &             ACCESS='SEQUENTIAL', STATUS='OLD')
	     READ(20,*) Mu00Data(I1,J), Mu10Data(I1,J), Mu01Data(I1,J),
     &		        Mu11Data(I1,J), Mu20Data(I1,J), Mu02Data(I1,J),
	&				Mu21Data(I1,J), Mu12Data(I1,J), Mu22Data(I1,J),
	&				Mu30Data(I1,J), Mu31Data(I1,J) 		 

	  ENDDO
	  CLOSE(UNIT=20)
 	ENDDO


      call FFSQP(nparam,nf,nineqn,nineq,neqn,neq,mode,iprint,
     *           miter,inform,bigbnd,eps,epseqn,udelta,bl,bu,x,f,g,
     *           iw,iwsize,w,nwsize,FCN,cntr,grobfd,grcnfd)


	OPEN(UNIT=20, FILE='param.dat', FORM='FORMATTED',

     &      ACCESS='SEQUENTIAL', STATUS='UNKNOWN')

      WRITE(20,10)x(1) 

	WRITE(20,10)x(2)

	WRITE(20,10)x(3)

	WRITE(20,10)x(4)

	WRITE(20,10)x(5)

	WRITE(20,10)x(6)


   10	FORMAT(E24.16)

      STOP
      END

*****************************************************************

      SUBROUTINE cntr(Nvar,jjj,x,gj)

*	cntr is the subroutine for the constraints.  The
*	constraints are listed in the following order:   
*	Nonlinear inequality constraints, linear inequality 
*	constraints, nonlinear equality constraints, and linear 
*	equality constraints.
*
*     input:   Nvar - number of parameters
*	         jjj  - indicates the jjj_th constraint
*	         x    - Nvar-dimensional vector of parameters
*	       
*     output:  gj   - the jjj_th constraint
*

      INTEGER Nvar,jjj
      REAL*8 x(*),gj


      RETURN
      END

******************************************************************
*
*       SUBROUTINE FCN(Nvar,j,THETA,Phi)
*
*       Nvar = number of parameters to be determined
*       THETA = vector of length N, parameters to be determined
*       Phi = objective function to be optimized
*
*       This subroutine simulates the operation of an industrial
*       crystallizer scaled-up from the experimental batch cooling
*       crystallizer described in S. M. Miller's Ph.D. thesis
*       published at the University of Texas at Austin in 1993.
*
*
*       Date:    February 12, 1998
*       Authors: Serena H. Chung and Richard D. Braatz
*                Department of Chemical Engineering
*                University of Illinois at Urbana-Champaign
*	  Modified:March 4, 2000

*	  By:      David L. Ma and Richard D. Braatz
*

      SUBROUTINE FCN(Nvar,j,THETA,Phi)

      INTEGER CSTEP, MSTEP, Nsets, Cinterval, Minterval
      PARAMETER (CSTEP=160,Nsets=1,MSTEP=6)
      PARAMETER (Cinterval=1,Minterval=30)

      INTEGER NN, NEQ,j,Nvar,MXPARM,Ndata
      PARAMETER (NN=CSTEP+1,NEQ=19,MXPARM=50,Ndata=Nsets)
      REAL*8 THETA(*), Phi
      INTEGER kfinal, I, Iteration
      INTEGER Norder, LDA, LDB, IPATH
      PARAMETER(Norder=3, LDA=3, LDB=3,IPATH=1)
      REAL*8 T, Y(NEQ)
      REAL*8 delt, tfinal
      REAL*8 mu00
      REAL*8 cell_length, Msolv, kv, ka, UA, densityc, densitys
      REAL*8 r0, alpha, g1, g2, kg1, kg2, b, kb
      
	REAL*8 moment00(NN), moment10(NN), moment01(NN)
	REAL*8 moment11(NN), moment20(NN), moment02(NN)
	REAL*8 moment21(NN), moment12(NN), moment22(NN)
      REAL*8 moment30(NN), moment31(NN)
      
	REAL*8 time (NN), concentration(NN)
	REAL*8 concentration_measured(NN)
      REAL*8 seed_moment10(NN), seed_moment01(NN)
	REAL*8 seed_moment11(NN), seed_moment20(NN)
	REAL*8 seed_moment02(NN), seed_moment21(NN)
	REAL*8 seed_moment12(NN)
      REAL*8 temperature(NN), relsatn(NN)
      REAL*8 Temp, Csat

      REAL*8 ConcData(Ndata,CSTEP), Mu00Data(Nsets,MSTEP)
      REAL*8 Mu10Data(Nsets,MSTEP), Mu01Data(Nsets,MSTEP)
      REAL*8 Mu11Data(Nsets,MSTEP), Mu20Data(Nsets,MSTEP)
	REAL*8 Mu02Data(Nsets,MSTEP), Mu12Data(Nsets,MSTEP)
      REAL*8 Mu21Data(Nsets,MSTEP), Mu22Data(Nsets,MSTEP)
      REAL*8 Mu30Data(Nsets,MSTEP), Mu31Data(Nsets,MSTEP)
	REAL*8 weight_conc, weight_mu00, weight_mu10
      REAL*8 weight_mu01, weight_mu11, weight_mu20
	REAL*8 weight_mu02, weight_mu21, weight_mu12
      REAL*8 weight_mu22, weight_mu30, weight_mu31

      REAL*8 AA(3,3), BB(3), gamma(3), lmin, lmax
      REAL*8 mass_seed, mean_seed, width

     
*lsodes' parameters
      
      INTEGER itol, iopt, itask, istate, mf
      INTEGER lrw, liw, iwork(200)
      REAL*8  rtol, atol, rwork(3800)      
            
      EXTERNAL MOMENTS, MOMENTSJ, DLSARG
      COMMON /GROWTH_DATA1/kg1, g1
	COMMON /GROWTH_DATA2/kg2, g2
      COMMON /BIRTH_DATA/kb, b 
      COMMON /EXP_DATA/r0, alpha, mu00, UA, Msolv
      COMMON Iteration
      COMMON /CONCEN/ConcData
      COMMON /MU00/Mu00Data
      COMMON /MU10/Mu10Data
      COMMON /MU01/Mu01Data
      COMMON /MU11/Mu11Data 
      COMMON /MU20/Mu20Data
	COMMON /MU02/Mu02Data
      COMMON /MU21/Mu21Data
      COMMON /MU12/Mu12Data 
      COMMON /MU22/Mu22Data
	COMMON /MU30/Mu30Data 
      COMMON /MU31/Mu31Data

      COMMON /WEIGHT/weight_conc,
     &			 weight_mu00, weight_mu10, weight_mu01,
     &			 weight_mu11, weight_mu20, weight_mu02,
     &			 weight_mu21, weight_mu12, weight_mu22,
     &		     weight_mu30, weight_mu31 
 
     
      Phi=0.0D0
      error_conc=0.0D0
      error_trans=0.0D0
 
      DO 777 Iteration=1, Ndata
	
	
*Simulation parameters
*     controller time step in minutes
      delt=1.0D0
*     total time step
      kfinal=NN
*     final time in minutes
      tfinal=DFLOAT(NN)*delt

*Parameters for experimental set-up
*     cell length for spectrophotometer in millimeter
*     This was modified from that in Miller's thesis because
*     his value (2.0) did not agree with his simulation results
      cell_length=1.77D0
*     mass of solvent in grams, converted from 2000 gallons
      Msolv=7.57D6
*     volume shape factor (Appendix C in Miller)
      kv=1.0D0
*     area shape factor (Appendix C in Miller)
      ka=6.0D0
*     heat transfer coefficient multiplied by surface area
*     in calorie/minute/degree C 
*     density of solvent in g/cm^3 (solvent is water)
      densitys=1.0D0
*     density of crystal in g/cm^3 (Appendix C in Miller)
      densityc=2.11D0
*     seed size at nucleation
       r0=0.0D0
*     crystal density*volume shape factor,
*     in gram crystal/micron^3/particle
*     (alpha*L^3=mass of particle)
      alpha=kv*densityc*(1.0D-4)**3

*     Total mass of seed crystals (grams)	
      mass_seed = 230.0D0
*     Mean size of the seed crystals (microns)