#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#define NMAX 101

/*****************************************************************************
 *
 * Calculates an array of lambda pairs (lam1,lam2) for polar values of theta
 * and nu (nu**2 is the 2nd moment about the origin).  Thus:
 *     radius = nu, nusq = mu**2 = mu**2 + sigma**2,
 *     mu = radius * cos(theta),
 *     sigma = radius * sin(theta).
 * theta varies from pi - del_th to del_th,
 * radius varies from del_r to 1 - del_r.
 *
 *
 *****************************************************************************
 *****************************************************************************
 *
 * Harold J. Price 
 * October 10, 2001
 *
 *****************************************************************************
 *****************************************************************************/

int
main(void)
{



    int numx; /* number of x points for x integration */
    int i,j,m,num_iter;
    int num_mu,num_sig,n_mu,n_sig;
    
    double xn,xnsq,delx,a,b,lam1,lam2,lam1p,lam2p;
    double g1_n,g2_n;
    double j11,j12,j21,j22,det;
    double k11,k12,k21,k22;
    double f0_n,f1_n,f2_n,f_lamn;
    double mu,sigma,musq,sigsq,nusq;
    double eps;
    double del_mu,del_sig;
    double lambda_1[NMAX][NMAX];
    double lambda_2[NMAX][NMAX];
    double norm[NMAX][NMAX];
    double pofmu[NMAX][NMAX];

    FILE *fp1,*fp2,*fp3;

    eps = 1.0e-14;
    numx = 1000;
    num_mu = NMAX - 1;
    num_sig = NMAX - 1;
    num_iter = 25;

    a = -1.0;
    b = 1.0;
    delx = (b - a)/(numx - 1);
    del_mu = (b-a)/(num_mu -1);
    del_sig = (1.0 - 0.0)/(num_sig - 1);


    fp1 = fopen("lambda_1.fil","w");
    fp2 = fopen("lambda_2.fil","w");
    fp3 = fopen("f0_n.fil","w");

/* begin mu - nu loop : nusq = musq + sigsq */

    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        mu =  a + (n_mu - 1) * del_mu;
        musq = mu * mu;

        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {      
            sigma = (n_sig - 1) * del_sig;
            sigsq = sigma * sigma;
            nusq = musq + sigsq;

            if((sqrt(nusq) >= 1) || (sigma < 0.01))
            {
                     lam1 = 0.0;
                     lam2 = 0.0;
                     f0_n = 0.0;
                     goto end_loop;
            }

/* begin while loop to find lambdas */

            g1_n = 1.0;
            g2_n = 1.0;
            lam1 = -1.0;
            lam2 = 1.0;

            m = 1;
            while(((fabs(g1_n) >= eps) || (fabs(g2_n) >= eps)) && (m <= num_iter))
            {

/* begin integration loop */

                xn = a;
                xnsq = xn * xn;
                f_lamn = exp(-lam1 * (xn - mu) - lam2 * (xnsq - nusq)) * delx;
                f0_n = 0.5 * f_lamn;
                f1_n = -(xn - mu) * f_lamn;
                f2_n = -(xnsq - nusq) * f_lamn;
                g1_n = 0.5 * f1_n;
                g2_n = 0.5 * f2_n;
                j11 = -(xn - mu) * g1_n;
                j12 = -(xn - mu) * g2_n;
                j21 = -(xnsq - nusq) * g1_n;
                j22 = -(xnsq - nusq) * g2_n;

                for(i = 2;i <= (numx - 1);i++)
                {
                    xn = a + (i - 1) * delx;
                    xnsq = xn * xn;
                    f_lamn = exp(-lam1 * (xn - mu) - lam2 * (xnsq - nusq)) * delx;
                    f0_n += f_lamn;
                    f1_n = -(xn - mu) * f_lamn;
                    f2_n = -(xnsq - nusq) * f_lamn;
                    g1_n += f1_n;
                    g2_n += f2_n;
                    j11 += -(xn - mu) * f1_n;
                    j12 += -(xn - mu) * f2_n;
                    j21 += -(xnsq - nusq) * f1_n;
                    j22 += -(xnsq - nusq) * f2_n;
/*            
            printf("i,f1_n,f2_n,xn: %d %e %e %e\n",i,f1_n,f2_n,xn);
*/
                }

                xn = b;
                xnsq = xn * xn;
                f_lamn = exp(-lam1 * (xn - mu) - lam2 * (xnsq - nusq)) * delx;
                f0_n += 0.5 * f_lamn;
                f1_n = -(xn - mu) * f_lamn;
                f2_n = -(xnsq - nusq) * f_lamn;
                g1_n += 0.5 * f1_n;
                g2_n += 0.5 * f2_n;
                j11 += 0.5 * -(xn - mu) * f1_n;
                j12 += 0.5 * -(xn - mu) * f2_n;
                j21 += 0.5 * -(xnsq - nusq) * f1_n;
                j22 += 0.5 * -(xnsq - nusq) * f2_n;

/* end integration loop; begin Newton-Raphson */

    /* k's are the elements of the invers Jacobian matrix with elements j */

                det = j11 * j22 - j12 * j21;
                k11 = j22/det;
                k12 = -j12/det;
                k21 = -j21/det;
                k22 = j11/det;

                lam1p = lam1 - (k11 * g1_n + k12 * g2_n);
                lam2p = lam2 - (k21 * g1_n + k22 * g2_n);

/* end Newton-Raphson */
/*
                printf("1,2,1p,2p: %e, %e, %e, %e \n",lam1,lam2,lam1p,lam2p);
                printf("m,g1_n,g2_n: %d %e %e\n",m,g1_n,g2_n);
                printf("m,f0_n, : %d %e\n\n",m,f0_n);
*/

                lam1 = lam1p;
                lam2 = lam2p;
         
                m += 1;
    
            }/* end of while loop */
/*
            printf("n_mu,n_sig,lam1,lam2 : %d %d %e %e \n",n_mu,n_sig,lam1,lam2);
*/

end_loop:

            fprintf(fp1,"%e %e %e\n",mu,sigma,lam1);
            fprintf(fp2,"%e %e %e\n",mu,sigma,lam2);
            if(f0_n != 0.0)
            {
                fprintf(fp3,"%e %e %e\n",mu,sigma,log(f0_n));
            }
            else
            {
                fprintf(fp3,"%e %e %e\n",mu,sigma,f0_n);
            }


            lambda_1[n_mu][n_sig] = lam1;
            lambda_2[n_mu][n_sig] = lam2;
            if(f0_n != 0,0) norm[n_mu][n_sig] = log(f0_n);


        }/* end n_sig loop */
    
        printf("n_mu,n_sig %d %d\n\n",n_mu,n_sig - 1);
       
    }/* end n_mu loop */ 

    fclose(fp1);
    fclose(fp2);
    fclose(fp3);

} /* end of main */