/******************************************************************************
*   This PL/1 subroutine calculates the ion densities in local chemical 
*   equilibrium for an arbitray number of constituents under consideration 
*   of ionization, recombination and charge exchange processes.
*   The program is based on the equilibrium condition for the production and 
*   loss rates for the ion species i
*   
*   qi +Σl(≠i)ΣmΣp ki,l,m,p.nl.Nm  = αi.ni.(ni +Σl(≠i)nl) + ni.ΣmΣp ki,m,p.Nm  
*   where
*
*   qi : primary ion production rate for species i
*   αi : recombination coefficient for ion i
*   Nm : neutral density of species m
*   ki,m,p: rate coefficent for charge exchange of ion i with neutral m
*   ki,l,m,p: rate coefficient for reaction of ion l and neutral m resulting in ion i
*
*   (for the charge exchange reactions, the index l sums over all ionized 
*   constituents,  the index m over all corresponding neutral constituents, 
*   and the index p sums over  possible multiple reaction channels (e.g. 
*   N++O2 → N+O2+  and N++O2 → O+NO+) ).
*
*   This is a quadratic equation for the ion density ni which can be 
*   readily solved to yield
*
*   ni = -Ai/2 +√[Ai2/4 +qi/αi +1/αi.Σl(≠i)ΣmΣp ki,l,m,p.nl.Nm ]  , 
*
*   where
*
*   Ai = Σl(≠i)nl + 1/αi.ΣmΣp ki,m,p.Nm .
*
*   As this solution is a function of the other ion densities nl, this 
*   constitutes a system of M (the number of constituents) equations 
*   that can be solved by iteration (which is what this program does).
*
*   It is recommended that this program uses internally an extended 
*   numerical accuracy, because of the possibility that near equal numbers 
*   are being subtracted from each other.
*
*
*  Input parameters are:
*    NC:  number of constituents
*    NMC: maximum number of charge exchange channels for any constituent
*    EPSI: desired relative accuracy
*    NEUDEN: neutral gas densities (in cm-3) (vector of length NC)
*    QP: ionization rates (in cm-3sec-1) (vector of length NC)
*             (note that the original vector content is modified in the program)
*    AL: recombination coefficient (in cm3sec-1) (vector of length NC)
*    KR: charge exchange coefficients (in cm3sec-1)(array of dimension NC*NC*NMC)
*    KRL: reaction channels for the production by charge exchange 
*       (=1: reaction is taken into account; otherwise =0) (array of dimension NC*NC*NC*NMC)
*
*  Output parameter is:
*    IONDEN: ion densities (in cm-3) (vector of length NC)
*
*    (other units can be used, but must obviously be consistent).
*
*   This program has worked well and correct for me in the past (using IBM PL/1 compilers),  
*   but I can't guarantee this for all compilers, especially if the code is translated
*   into a different language.
*
*   If there are any problems, please let me know and I mention the issue here.
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*                                               Thomas Smid, April 2007
*********************************************************************************/

 IONDENS:PROCEDURE(NC,NMC,EPSI,NEUDEN,QP,AL,KR,KRL,IONDEN) ;
 DEFAULT RANGE(A:Z) FLOAT DECIMAL VALUE (DEC FLOAT(33)) ;
 DCL (DEN(*),L(*),NREL(*),QN(*),NDEN(*)) CTL DEC FLOAT(33) ;
 DCL (NQN(*),ID(*)) CTL DEC FIXED(7,0) ;
 DCL (NEUDEN(*),QP(*),AL(*),KR(*,*,*)) CTL DEC FLOAT(16) ;
 DCL KRL(*,*,*,*) CTL DEC FIXED(3,0) ;
 DCL (I,J,K,M,N,NC,DI,NMC,SNQN,NUMIT) DEC FIXED(7,0) ;
 DCL IONDEN(*) CTL DEC FLOAT(16) ;
 ON UNDERFLOW ; 


 ALLOCATE L(NC),NREL(NC),NQN(NC),DEN(NC),NDEN(NC) ;
 DEN=0; IONDEN=0; NDEN=0;
 LLIM=1E-33/EPSI ;   /* 1E-33 should be changed to the available precision if different */


 /* Determination of total number of charge exchange gain reactions */
 NQN=0 ;  
 DO M=1 TO NC ;    
 DO I=1 TO NC ; 
 DO J=1 TO NC ; 
 DO K=1 TO NMC ;
 NQN(M)=NQN(M)+KRL(M,I,J,K) ;  
 END ;
 END ;
 END ;
 END ;
 SNQN=SUM(NQN) ;
 ALLOCATE QN(SNQN),ID(SNQN) ;


 /* Calculation of charge exchange loss rates and initial ion densities */
 L=0 ;
 DO I=1 TO NC ; 
 DO J=1 TO NC ; 
 DO K=1 TO NMC ;
 L(I)=L(I)+KR(I,J,K)*NEUDEN(J) ;
 END ;
 END ;
 L(I)=L(I)/2/AL(I)     ;
 QP(I)=QP(I)/AL(I) ;
 DEN(I)=-L(I)+SQRT(QP(I)+L(I)**2);
 IF L(I)>0 THEN DO ;
 QPS=QP(I)/L(I)**2 ;
 IF QPS<LLIM|DEN(I)<LLIM THEN CALL ERPRINT ;
 END ;
 END ;

  /* Calculation of charge exchange gain rates */
 N=1; QN=0; ID=0;
 DO M=1 TO NC ; 
 DO I=1 TO NC ; 
 DO J=1 TO NC ; 
 DO K=1 TO NMC ;
 IF KRL(M,I,J,K)=1 THEN DO ;
 QN(N)=KR(I,J,K)*NEUDEN(J)/AL(M);
 ID(N)=I ;
 N=N+1 ;
 END ;
 END ;
 END ;
 END ;
 END ;

 NDEN=DEN; NUMIT=0;


 SDEV=1; JJ=0; DI=1;
 ITER: DO WHILE(SDEV>=EPSI) ;     /* Iteration loop */
 NUMIT=NUMIT+1; JJ=JJ+1; NDEN=0;
 ELDEN=SUM(DEN) ;
 N=1;
 DO I=1 TO NC ; 
 DO J=1 TO NQN(I) ;
 NDEN(I)=NDEN(I)+QN(N)*DEN(ID(N)) ;
 N=N+1 ;
 END ;
 SQA=-(ELDEN-DEN(I))/2-L(I) ;
 QPN=QP(I)+NDEN(I) ;
 SQAN=SQRT(QPN+SQA**2) ;
 NDEN(I)=SQA+SQAN ;
 IF SQA<0 THEN DO ;
 QPS=QPN/SQA**2 ;
 IF QPS<LLIM|NDEN(I)<LLIM THEN CALL ERPRINT ;
 END ;
 END ;


 DO I=1 TO NC ; 
 IF DEN(I)=0 THEN GO TO CONT1 ; 
 NREL(I)=NDEN(I)/DEN(I) ;
 CONT1:END ;


 SDEV=ABS(NREL(1)-1) ;
 DO I=2 TO NC ; 
 SDEV=MAX(SDEV,ABS(NREL(I)-1)) ;  /* check for convergence */
 END ;
 SDEV=SDEV*JJ ; 
 DEN=NDEN ;

 END ITER ;     /* End 'Iteration loop' */



 IF DI=1 THEN DO ;
 DI=-1; SDEV=1; 
 GO TO ITER ;    /* Additional test iteration after convergence is achieved */
 END ;

 IONDEN=NDEN ;

 FREE DEN,NDEN,L,QN,NREL,ID,NQN ;



 ERPRINT: PROCEDURE ;
 PUT EDIT('DESIRED ACCURACY CAN NOT BE REACHED') (A(35)) ;
 IONDEN=1E50 ;
 GO TO ENDE ;
 END ERPRINT ;



 ENDE:END IONDENS ;