#include "peg_pic.h"


/** Cross section of ion-neutral particle charge-exchange collisions **/

float sigma5(int ii, int k, float energy) {
	int	i, j, p;
	float 	deng;
	static float 	***sicx1, ***sicx2, ***sicx3;
	static int 	init_flag=1;
	char  	tchar[80];


	/**********  Initialization  **************/
	if(init_flag) {
		p=0;
		
		for (i=0; i<Nnum; i++) {
			for (j=0; j<Nnum; j++) {
				if (HH[i][j]) {
					p +=1;
				}
			}
		}

		if (p) {                /* allocate space for exact cross-sections */
			sicx1= (float ***)malloc(Nnum*sizeof(float **));
			sicx2= (float ***)malloc(Nnum*sizeof(float **));
			sicx3= (float ***)malloc(Nnum*sizeof(float **));

			for (i=0; i<Nnum; i++) {
				sicx1[i]= (float **)malloc(Nnum*sizeof(float *));
				sicx2[i]= (float **)malloc(Nnum*sizeof(float *));
				sicx3[i]= (float **)malloc(Nnum*sizeof(float *));

				for (j=0; j<Nnum; j++) {
					sicx1[i][j]= (float *)malloc(1000*sizeof(float));
					sicx2[i][j]= (float *)malloc(200*sizeof(float));
					sicx3[i][j]= (float *)malloc(NEMAX*sizeof(float));
				}
			}
		}   /* end of if (p) loop */


		for (i=0; i<Nnum; i++) {
			for (j=0; j<Nnum; j++) {
				if (HH[i][j]) {                 /* exact cross-sections */
					dec1(tchar, &p, 1, HH[i][j], 1, 0, 0);

					if (!p) {                    /* reading from the table */
						tcrossec(sicx1[i][j], tchar, 1000, 1.e-4);  // from 0.0 to 0.1
						tcrossec(sicx2[i][j], tchar, 200, 0.1);      // from 0.1 to 20
						tcrossec(sicx3[i][j], tchar, NEMAX, 1.);     // from 20 to NEMAX
					} else {                        /* calculating from the fit function */
						for (p=0; p<1000; p++) {
							deng= 1.e-4*p;
							sicx1[i][j][p] = (*fcrs)(deng);
						}

						for (p=0; p<200; p++) {
							deng= 0.1*p;
							sicx2[i][j][p] = (*fcrs)(deng);
						}

						for (p=0; p<NEMAX; p++) {
							sicx3[i][j][p] = (*fcrs)((float)p);
						}
					}

					if (dshort>0.) {
						for (i=0; i<Nnum; i++) {
							for (j=0; j<Nnum; j++) {
								for (p=0; p<1000; p++) {
									sicx1[i][j][p] *= dshort;
								}
								
								for (p=0; p<200; p++) {
									sicx2[i][j][p] *= dshort;
								}
								
								for (p=0; p<NEMAX; p++) {
									sicx3[i][j][p] *= dshort;
								}
							}
						}
					}
				} else if (dshort>0.) {              /* approximated cross-sections */
					achrgx[i][j] *=dshort;
					bchrgx[i][j] *=dshort;
				}
			}
		}      /* end of for (i=0; i<Nnum; i++) loop */

		init_flag =0;

	}  /* end of initialization */


	if (HH[ii][k]) {
		if (energy<0.1) {
			i= 10000*energy;
			return(sicx1[ii][k][i]);
		} else if (energy<20) {
			i= 10*energy;
			return(sicx2[ii][k][i]);
		} else if (energy<NEMAX) {
			i= energy;
			return(sicx3[ii][k][i]);
		} else {
			return(0.);
		}
	} else {
		return(achrgx[ii][k] + bchrgx[ii][k]/(sqrt(energy) +1e-30));
	}
}
/* END of sigma5 */

/************************************************/
