/*

RSII soln 5d

*/

/*
 
 Direct solver; compute Cheb internally; variable precision
 
*/



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

/*
 
 REAL       float   double  long double
 REALTYPE   0       1       2
 
 */

#define REAL      double 
#define REALTYPE 1


#define LAPACK  1 


#define Power(x,n) powl(x,n)
#define Sqrt(x)    sqrtl(x)

#define PI 3.1415926535897932384626433832795028841971693993751


#define LICHCUT 0.


#define rhval 0.5

REAL rh = rhval ;

#define coeffAterm  +1.



#define Lambda -4.


#define EvenR 0
#define EvenA 0


#define NMAX  100
#define NMAX2 40000
#define NMAX3 50000000


REAL multval = -1. ;


int    NR,NA ;

REAL   patch[NMAX][NMAX][2] ;
long   mask[5][NMAX][NMAX] ;
int    vecpos[NMAX2][3] ;
long   vecsize, sizeLich ;

REAL   vecR[NMAX2], vecG[NMAX2], vecbuf[NMAX2], vecbuf2[NMAX2], vecRstore[NMAX2] ;

int    LichRow[NMAX3], LichCol[NMAX3] ;
REAL   Lich[NMAX3] ;

REAL   *Lmat ;



REAL DRmat[NMAX][NMAX], DR2mat[NMAX][NMAX], DAmat[NMAX][NMAX], DA2mat[NMAX][NMAX], DA3mat[NMAX][NMAX] ;

REAL gmet[5][NMAX][NMAX] ; 
REAL RF[5][NMAX][NMAX] ;
REAL phi[NMAX][NMAX], Ricci[NMAX][NMAX], R2sqr[NMAX][NMAX], R4sqr[NMAX][NMAX] ;

REAL dvr[5][NMAX][NMAX], dva[5][NMAX][NMAX], dvrr[5][NMAX][NMAX], dvaa[5][NMAX][NMAX], dvaaa[5][NMAX][NMAX], dvra[5][NMAX][NMAX] ;



void initpatches(void), init(void), save(int) ;
void calcRicciFlow(REAL *), calcLich(void) ;

REAL calcPhi(void) ;

void getG(void), storeG(void), updateG(REAL, REAL *) ;


REAL valarr[5] ;
FILE *fileptr, *fileptrphi, *fileptrRicci, *fileptrR2sqr, *fileptrR4sqr, *fileptrDerivs ;
REAL DT ;


void unpackLich(void) ;

void linsolve(REAL *,REAL *), linsolvechk(REAL *,REAL *) ;
REAL norm(REAL *) ;


void setupCheb(void) ;

int status ;



int main (int argc, const char * argv[]) 
{
	int iter ;
	REAL rad,multiplier,residual,Rerr,phimax ;
	clock_t firsttime, starttime, time1, time2 ;  
    long Lsize, posii, posjj ;

    
#if REALTYPE == 2
    
    printf("\n\nLONG DOUBLE\n\n") ;
    
#elif REALTYPE == 1    
    
    printf("\n\nDOUBLE\n\n") ;
    
#else    
    
    printf("\n\nFLOAT\n\n") ;
    
#endif 
    
    
    
    if(LAPACK!=0 && REALTYPE!=1) {
        
        printf("\n\nError - LAPACK must use *double*\n\n") ; exit(1) ;
        
    }
    
	
	fileptr = fopen("data.out","w") ; fileptrphi = fopen("phi.out","w") ; 
    fileptrRicci = fopen("Ricci.out","w") ;
    fileptrR2sqr = fopen("R2sqr.out","w") ;
    fileptrR4sqr = fopen("R4sqr.out","w") ;
    fileptrDerivs = fopen("derivs.out","w") ;
	
	initpatches() ;
		
	init() ; 
    
    
    
	printf("\n\n ** Initial values:  NR = %d  NA = %d   rh = %f    Lambda = %f \n\n",NR,NA,rhval,Lambda) ;	
	
    
    setupCheb() ;
    

    printf("\n\nr coords:\n\n") ;
#if REALTYPE == 2                                         
    for(posii=0;posii<NR;posii++) printf("   %Lf\n",patch[posii][0][0]) ;   
#else                                        
    for(posii=0;posii<NR;posii++) printf("   %f\n",patch[posii][0][0]) ;                               
#endif     

    printf("\n\na coords:\n\n") ;
#if REALTYPE == 2                                         
    for(posii=0;posii<NA;posii++) printf("   %Lf\n",patch[0][posii][1]) ;   
#else                                        
    for(posii=0;posii<NA;posii++) printf("   %f\n",patch[0][posii][1]) ;                               
#endif     

    printf("\n\n") ;
    
    
    Lsize = vecsize*(vecsize+1)*sizeof(REAL) ;
    Lmat  = malloc(Lsize) ;
    
    if(Lmat==NULL) {
    
        printf("\n\nCouldn't allocate memory for Lich op matrix\n\n") ;
    
    } else {
        
        printf("\n\nAllocated Lich memory:  %ld  bytes\n\n     %ld  Mb\n\n",Lsize,Lsize/(1000000)) ;
        
    }
    
	printf("\n\nNewton solver:\n\n") ;
	
	multiplier = multval ;

#if REALTYPE == 2
 
            printf("multiplier %Lf\n",multiplier) ;
            
#else    
    
            printf("multiplier %f\n",multiplier) ;
            
#endif 
            
            
	calcRicciFlow(vecR) ;
	
	firsttime = clock() ;
    
    status = 1 ;
	
	for(iter=0;iter<20;iter++) {
		
#if REALTYPE == 2
        
        printf("Iteration = %d ;  multiplier = %Lf  status = %d   rh = %Lf\n",iter,multiplier,status,rh) ;
        
#else    
        
        printf("Iteration = %d ;  multiplier = %f  status = %d    rh = %f\n",iter,multiplier,status,rh) ;
        
#endif 

				
		starttime = clock() ;
				
		printf("  <- ") ; fflush(stdout) ;

		calcRicciFlow(vecR) ;
	
		calcLich() ;	
        
        unpackLich() ;
		
		printf("  ->\n") ;   fflush(stdout) ;
		
		
		time1 = clock() ;
		
		linsolve(vecbuf,vecR) ;
		
               
        unpackLich() ;
		linsolvechk(vecbuf,vecR) ;


    
        for(posii=0;posii<vecsize;posii++) vecRstore[posii] = vecR[posii] ;

        

		time2 = clock() ;
 		
		
		getG() ;
		
		updateG(multiplier,vecbuf) ;
		
		storeG() ;		
		
		residual = norm(vecbuf) ;
		
		calcRicciFlow(vecR) ;

		Rerr = norm(vecR) ;		
        
        phimax = calcPhi() ;

        save(0) ;	

        
#if REALTYPE == 2
        
        printf("   Update = %10.3Le  Rerr = %10.3Le  Phimax = %10.3Le  -- Timings (sec; Lich, Solve): %10.3f %10.3f\n",residual,Rerr,phimax,(1.*(time1 - starttime))/(1.*CLOCKS_PER_SEC),(1.*(time2 - time1))/(1.*CLOCKS_PER_SEC)) ;
        
#else    
        
        printf("   Update = %10.3e  Rerr = %10.3e  Phimax = %10.3e  -- Timings (sec; Lich, Solve): %10.3f %10.3f\n",residual,Rerr,phimax,(1.*(time1 - starttime))/(1.*CLOCKS_PER_SEC),(1.*(time2 - time1))/(1.*CLOCKS_PER_SEC)) ;
        
#endif 

        
        if(residual<1e-4 && status==0) {
			status = 1 ;
		}
         
		
		if(residual<1e-9) {
			printf("*** Finished\n") ; break ;
		}
				
	}
		
	printf("   Total Time: %10.3f (min)\n",(1./60.*(clock() - firsttime))/(1.*CLOCKS_PER_SEC)) ;
	
    
    unpackLich() ;
	save(1) ;	
	
    
    
	free(Lmat) ;
	
	fclose(fileptr) ;  fclose(fileptrphi) ; 
    
    fclose(fileptrRicci) ; fclose(fileptrR2sqr) ; fclose(fileptrR4sqr) ; 
    
    fclose(fileptrDerivs) ;

	
	return 0;
	
}



	

void calcRicciFlow(REAL *vecRptr)
{
	int ii, jj, var, kk ;
	long pos ;
	
	REAL val ;
	
	REAL EqX[5], EqBrn[5] ;
	
	REAL r,a ;
	
	REAL T, Tr, Ta, Trr, Taa, Tra ;
	REAL A, Ar, Aa, Arr, Aaa, Ara ;
	REAL B, Br, Ba, Brr, Baa, Bra ;
	REAL F, Fr, Fa, Frr, Faa, Fra ;
	REAL S, Sr, Sa, Srr, Saa, Sra ;
	
	REAL L, Lr, La, Laa, Lrr, Lra, f, fr, frr, g, ga, gaa, r2, a2 ;
	
    
	
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			for(var=0;var<5;var++) {
				
				dvr[var][ii][jj] = 0. ;
				for(kk=0;kk<NR;kk++) dvr[var][ii][jj] += DRmat[ii][kk]*gmet[var][kk][jj] ;
				
				dva[var][ii][jj] = 0. ;
				for(kk=0;kk<NA;kk++) dva[var][ii][jj] += DAmat[jj][kk]*gmet[var][ii][kk] ;
				
				dvrr[var][ii][jj] = 0. ;
				for(kk=0;kk<NR;kk++) dvrr[var][ii][jj] += DR2mat[ii][kk]*gmet[var][kk][jj] ;
				
				dvaa[var][ii][jj] = 0. ;
				for(kk=0;kk<NA;kk++) dvaa[var][ii][jj] += DA2mat[jj][kk]*gmet[var][ii][kk] ;
                
                dvaaa[var][ii][jj] = 0. ;
				for(kk=0;kk<NA;kk++) dvaaa[var][ii][jj] += DA3mat[jj][kk]*gmet[var][ii][kk] ;
				
			}
			
		}	
		
	}
    
    
    
 
	
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			for(var=0;var<5;var++) {
				
				dvra[var][ii][jj] = 0. ;
				for(kk=0;kk<NR;kk++) dvra[var][ii][jj] += DRmat[ii][kk]*dva[var][kk][jj] ;
				
			}
			
		}	
		
	}
	
	
	
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			
			r=patch[ii][jj][0] ; a=patch[ii][jj][1] ;
			
			
			
			L = Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,-2) ;
			
			Lr = (-4*r*rh)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,3);
			
			La = (-4*a)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,3);
			
			Lrr = (4*rh*(1 - Power(a,2) + rh + 5*Power(r,2)*rh))/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			Laa = (4*(1 + 5*Power(a,2) + rh - Power(r,2)*rh))/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			Lra = (24*a*r*rh)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			
			
			f = 1 - Power(r,2) ;
			
			fr = -2*r ;
			
			frr = -2 ;
			
			g = 2. - a*a ;
			
			ga = -2.*a ;
			
			gaa = -2. ;
			
			r2 = r*r ;
			
			a2 = a*a ;
			
			
			T = gmet[0][ii][jj] ; Tr = dvr[0][ii][jj] ; Ta = dva[0][ii][jj] ;
			A = gmet[1][ii][jj] ; Ar = dvr[1][ii][jj] ; Aa = dva[1][ii][jj] ; 
			B = gmet[2][ii][jj] ; Br = dvr[2][ii][jj] ; Ba = dva[2][ii][jj] ; 
			F = gmet[3][ii][jj] ; Fr = dvr[3][ii][jj] ; Fa = dva[3][ii][jj] ; 
			S = gmet[4][ii][jj] ; Sr = dvr[4][ii][jj] ; Sa = dva[4][ii][jj] ; 	
			
			Trr = dvrr[0][ii][jj] ; Taa = dvaa[0][ii][jj] ; Tra = dvra[0][ii][jj] ;
			Arr = dvrr[1][ii][jj] ; Aaa = dvaa[1][ii][jj] ; Ara = dvra[1][ii][jj] ;
			Brr = dvrr[2][ii][jj] ; Baa = dvaa[2][ii][jj] ; Bra = dvra[2][ii][jj] ;
			Frr = dvrr[3][ii][jj] ; Faa = dvaa[3][ii][jj] ; Fra = dvra[3][ii][jj] ;
			Srr = dvrr[4][ii][jj] ; Saa = dvaa[4][ii][jj] ; Sra = dvra[4][ii][jj] ;			
            
            
            
			
			if(ii==0 && EvenR != 1) {
				
				EqX[0] = Tr ; EqX[1] = Ar ; EqX[2] = Br ; EqX[3] = Fr ; EqX[4] = Sr ;
				
			} else if(jj==0 && EvenA != 1) {
				
				EqX[0] = Ta ; EqX[1] = Aa ; EqX[2] = Ba ; EqX[3] = Fa ; EqX[4] = Sa ;
				
			} else if(jj==NA-1) {
			
#include "Brn.c"				
                
                 
			} else {
			
#include "Eqn.c"
                
        
                if(ii==NR-2) {
                
                    EqX[1] += coeffAterm*dvra[1][NR-1][jj] ;
                    
                }
                
               
				
			}
			
			for(var=0;var<5;var++) {
				
				val = EqX[var] ;
					
				if( (pos=mask[var][ii][jj]) != -1) {
					
					vecRptr[pos] = val ;
					
					RF[var][ii][jj] = val ;
					
				}	
				
			}
			
		}
	}

    
	
}




void calcLich(void)
{
	int ii, jj, ii2, jj2, var, var2 ;
	
	long posii, posjj, posLich ;
	
	REAL val ;
	
	REAL coeffrr, coeffaa, coeffra, dEqXr[5][5], dEqXa[5][5], dEqT[5], dEqA[5], dEqB[5], dEqF[5], dEqS[5] ;
	
	REAL dBrnT[5], dBrnA[5], dBrnB[5], dBrnF[5], dBrnS[5] ;
	
	REAL r,a ;
	
	REAL T, Tr, Ta, Trr, Taa, Tra ;
	REAL A, Ar, Aa, Arr, Aaa, Ara ;
	REAL B, Br, Ba, Brr, Baa, Bra ;
	REAL F, Fr, Fa, Frr, Faa, Fra ;
	REAL S, Sr, Sa, Srr, Saa, Sra ;
	
	REAL L, Lr, La, Laa, Lrr, Lra, f, fr, frr, g, ga, gaa, r2, a2 ;
    
    
	posLich = 0 ;
	
	
	for(posii=0;posii<vecsize;posii++) {
		
		var = vecpos[posii][0] ;
		ii = vecpos[posii][1] ;
		jj = vecpos[posii][2] ;
		
		if(ii==0 && EvenR != 1) {
			
			for(posjj=0;posjj<vecsize;posjj++) {
				
				var2 = vecpos[posjj][0] ;
				ii2 = vecpos[posjj][1] ;
				jj2 = vecpos[posjj][2] ;
				
				if(var==var2 && jj==jj2) {	
					
					val = DRmat[0][ii2] ;
					
					LichRow[posLich] = posii ;
					LichCol[posLich] = posjj ;
					Lich[posLich]    = val ;
					
					posLich += 1 ;
					
					if(posLich>=NMAX3) {
						printf("posLich too large!\n") ; exit(1) ;
					}				
					
				}
				
			}
			
			
		} else if(jj==0 && EvenA != 1) {
			
			for(posjj=0;posjj<vecsize;posjj++) {
				
				var2 = vecpos[posjj][0] ;
				ii2 = vecpos[posjj][1] ;
				jj2 = vecpos[posjj][2] ;
				
				if(var==var2 && ii==ii2) {	
					
					val = DAmat[0][jj2] ;
					
					LichRow[posLich] = posii ;
					LichCol[posLich] = posjj ;
					Lich[posLich]    = val ;
					
					posLich += 1 ;
					
					if(posLich>=NMAX3) {
						printf("posLich too large!\n") ; exit(1) ;
					}				
					
				}
				
			}
			
		}  else {
			
			
			r=patch[ii][jj][0] ; a=patch[ii][jj][1] ;
			
			
			L = Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,-2) ;
			
			Lr = (-4*r*rh)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,3);
			
			La = (-4*a)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,3);
			
			Lrr = (4*rh*(1 - Power(a,2) + rh + 5*Power(r,2)*rh))/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			Laa = (4*(1 + 5*Power(a,2) + rh - Power(r,2)*rh))/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			Lra = (24*a*r*rh)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			
			f = 1 - Power(r,2) ;
			
			fr = -2*r ;
			
			frr = -2 ;
			
			g = 2. - a*a ;
			
			ga = -2.*a ;
			
			gaa = -2. ;
			
			r2 = r*r ;
			
			a2 = a*a ;
			
			
			T = gmet[0][ii][jj] ; Tr = dvr[0][ii][jj] ; Ta = dva[0][ii][jj] ;
			A = gmet[1][ii][jj] ; Ar = dvr[1][ii][jj] ; Aa = dva[1][ii][jj] ; 
			B = gmet[2][ii][jj] ; Br = dvr[2][ii][jj] ; Ba = dva[2][ii][jj] ; 
			F = gmet[3][ii][jj] ; Fr = dvr[3][ii][jj] ; Fa = dva[3][ii][jj] ; 
			S = gmet[4][ii][jj] ; Sr = dvr[4][ii][jj] ; Sa = dva[4][ii][jj] ; 	
			
			Trr = dvrr[0][ii][jj] ; Taa = dvaa[0][ii][jj] ; Tra = dvra[0][ii][jj] ;
			Arr = dvrr[1][ii][jj] ; Aaa = dvaa[1][ii][jj] ; Ara = dvra[1][ii][jj] ;
			Brr = dvrr[2][ii][jj] ; Baa = dvaa[2][ii][jj] ; Bra = dvra[2][ii][jj] ;
			Frr = dvrr[3][ii][jj] ; Faa = dvaa[3][ii][jj] ; Fra = dvra[3][ii][jj] ;
			Srr = dvrr[4][ii][jj] ; Saa = dvaa[4][ii][jj] ; Sra = dvra[4][ii][jj] ;			
            
            
            if(jj!=NA-1) {
			
#include "Licheq.c"
                
            } else {
                
#include "Lichbrn.c"       
                              
            }
			
			for(posjj=0;posjj<vecsize;posjj++) {
				
				var2 = vecpos[posjj][0] ;
				ii2 = vecpos[posjj][1] ;
				jj2 = vecpos[posjj][2] ;
				
				val = 0. ;
				
				if( ii==ii2 && jj==jj2 ) {
					
					switch(var2) {
							
						case(0) :
							
							val += dEqT[var] ;
							
							break ;
							
						case(1) :
							
							val += dEqA[var] ;
							
							break ;
							
						case(2) :
							
							val += dEqB[var] ;
							
							break ;
							
						case(3) :
							
							val += dEqF[var] ;
							
							break ;
							
						case(4) :
							
							val += dEqS[var] ;
							
							break ;
							
					}
					
				}
				
				
				if( jj==jj2 ) {
					
					if(var==var2)  val += coeffrr*DR2mat[ii][ii2] ;
					
					val += dEqXr[var][var2]*DRmat[ii][ii2] ;
					
				}
				
				
				
				if( ii==ii2 ) {
					
					if(var==var2)  val += coeffaa*DA2mat[jj][jj2] ;
					
					val += dEqXa[var][var2]*DAmat[jj][jj2] ;
					
				}
				
				if(var==var2)  {
                    
/* WARNING!! Modified here - no need to maintain sparseness */
											
						val += coeffra*DRmat[ii][ii2]*DAmat[jj][jj2] ;
                                     
					
				}
                
                
    /* MODIFIED - Arz and Axis  bc! */
                
                if(ii==NR-2 && jj<NA-1 && var==1 && var2==1) {
                    
                    val += coeffAterm*DRmat[NR-1][ii2]*DAmat[jj][jj2] ;
                    
                }
             
                
               /*  End modification */
                
                
				
				if(fabs(val)>LICHCUT) {	
					
					LichRow[posLich] = posii ;
					LichCol[posLich] = posjj ;
					Lich[posLich]    = val ;
					
					posLich += 1 ;
					
					if(posLich>=NMAX3) {
						printf("posLich too large!\n") ; exit(1) ;
					}				
					
					
				}
				
				
			}
			
			
		}
		
	}
	
	
	sizeLich = posLich ;
	
	
	printf("\nAnalytic: Size of Lich:  sizeLich = %ld\n",sizeLich) ;
	
}




REAL calcPhi(void)
{
	int ii, jj, var, kk ;
	
    REAL phimax ;
    
	REAL r,a ;
    
	REAL T, Tr, Ta, Trr, Taa, Tra ;
	REAL A, Ar, Aa, Arr, Aaa, Ara ;
	REAL B, Br, Ba, Brr, Baa, Bra ;
	REAL F, Fr, Fa, Frr, Faa, Fra ;
	REAL S, Sr, Sa, Srr, Saa, Sra ;
    
    REAL theta ;
	
	REAL hr, ha ;
	
	REAL L, Lr, La, Lrr, Laa, Lra, f, fr, frr, g, ga, gaa, r2, a2 ;
    
    REAL phival, Riccival, R2sqrval, R4sqrval ;
    
    REAL gd[5][5], gu[5][5], R2[5][5], C4[5][5][5][5] ;
    
    int i1,i2,i3,i4, s1, s2, s3, s4 ;
    
    phimax = 0. ;
    
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
            
			r=patch[ii][jj][0] ; a=patch[ii][jj][1] ;
            
			
			
			L = Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,-2) ;
			
			Lr = (-4*r*rh)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,3);
			
			La = (-4*a)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,3);
			
			Lrr = (4*rh*(1 - Power(a,2) + rh + 5*Power(r,2)*rh))/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			Laa = (4*(1 + 5*Power(a,2) + rh - Power(r,2)*rh))/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			Lra = (24*a*r*rh)/Power(-1 + Power(a,2) + (-1 + Power(r,2))*rh,4);
			
			
			f = 1 - Power(r,2) ;
			
			fr = -2*r ;
			
			frr = -2 ;
			
			g = 2. - a*a ;
			
			ga = -2.*a ;
			
			gaa = -2. ;
			
			r2 = r*r ;
			
			a2 = a*a ;
			
			
			T = gmet[0][ii][jj] ; Tr = dvr[0][ii][jj] ; Ta = dva[0][ii][jj] ;
			A = gmet[1][ii][jj] ; Ar = dvr[1][ii][jj] ; Aa = dva[1][ii][jj] ; 
			B = gmet[2][ii][jj] ; Br = dvr[2][ii][jj] ; Ba = dva[2][ii][jj] ; 
			F = gmet[3][ii][jj] ; Fr = dvr[3][ii][jj] ; Fa = dva[3][ii][jj] ; 
			S = gmet[4][ii][jj] ; Sr = dvr[4][ii][jj] ; Sa = dva[4][ii][jj] ; 	
            
            Trr = dvrr[0][ii][jj] ; Taa = dvaa[0][ii][jj] ; Tra = dvra[0][ii][jj] ;
			Arr = dvrr[1][ii][jj] ; Aaa = dvaa[1][ii][jj] ; Ara = dvra[1][ii][jj] ;
			Brr = dvrr[2][ii][jj] ; Baa = dvaa[2][ii][jj] ; Bra = dvra[2][ii][jj] ;
			Frr = dvrr[3][ii][jj] ; Faa = dvaa[3][ii][jj] ; Fra = dvra[3][ii][jj] ;
			Srr = dvrr[4][ii][jj] ; Saa = dvaa[4][ii][jj] ; Sra = dvra[4][ii][jj] ;			
            
            hr = (expl(-S - T)*(-4*f*expl(T)*(Lr*r*(-(Power(a,2)*Power(F,2)*g*Power(r,2)*
                                                      (-5*expl(S) + expl(A*Power(r,2) + S) + 2*expl(A*Power(r,2) + T))) + 
                                                    16*(-4*expl(Power(a,2)*B + S) + expl(Power(a,2)*B + A*Power(r,2) + S) + expl(A*Power(r,2) + T) + 
                                                        2*expl(Power(a,2)*B + A*Power(r,2) + T))*expl(A*Power(r,2) + S + T)) + 
                                              L*expl(S)*(Power(a,2)*Power(F,2)*g*Power(r,2)*(2 + 2*r*Sr + r*Tr - 2*expl(A*Power(r,2))) + 
                                                         32*A*Power(r,2)*expl(Power(a,2)*B + A*Power(r,2) + S + T) + 
                                                         16*(-2 - Power(a,2)*Br*r + Ar*Power(r,3) - 3*r*Sr + 2*expl(A*Power(r,2)))*expl(Power(a,2)*B + A*Power(r,2) + S + T)))
                                + r*(Power(a,2)*Power(F,3)*g*Power(r,3)*(2*a*ga*L*expl(T) + g*(4*L*expl(T) + a*La*(expl(S) + 2*expl(T)))) + 
                                     8*Power(a,2)*Power(F,2)*fr*g*L*Power(r,2)*(-1 + expl(A*Power(r,2)))*expl(S + T) + 
                                     16*F*r*(-(a*ga*L*(1 + 2*expl(Power(a,2)*B))*expl(A*Power(r,2) + T)) - 
                                             g*(2*L*(1 - a*Aa*Power(r,2) - a*Ta + 2*expl(Power(a,2)*B))*expl(A*Power(r,2) + T) + 
                                                a*La*(expl(Power(a,2)*B + S) + expl(Power(a,2)*B + A*Power(r,2) + S) + expl(A*Power(r,2) + T) + 
                                                      2*expl(Power(a,2)*B + A*Power(r,2) + T))))*expl(S + T) - 
                                     32*L*(a*Fa*g*r + 4*fr*(-1 + expl(A*Power(r,2)))*expl(Power(a,2)*B + S))*expl(A*Power(r,2) + S + 2*T))))/
            (8.*f*L*r*(-(Power(a,2)*Power(F,2)*g*Power(r,2)) + 16*expl(Power(a,2)*B + A*Power(r,2) + S + T))) ;
            
            ha = (expl(-S - T)*(Power(a,2)*Power(F,2)*Power(g,2)*Power(r,2)*
                                (2*L*expl(S)*(Power(a,2)*f*F + Power(a,2)*F*fr*r + 2*(-4 - 2*a*Sa - a*Ta + 4*expl(Power(a,2)*B))*expl(T)) + 
                                 a*(a*f*F*Lr*r*(expl(S) + 2*expl(T)) + 4*La*expl(S)*(expl(Power(a,2)*B + S) - 5*expl(T) + 2*expl(Power(a,2)*B + T)))) - 
                                8*g*expl(S + T)*(2*a*(4*La*expl(Power(a,2)*B + S)*
                                                      (expl(Power(a,2)*B + S) + expl(Power(a,2)*B + A*Power(r,2) + S) - 4*expl(A*Power(r,2) + T) + 
                                                       2*expl(Power(a,2)*B + A*Power(r,2) + T)) + 
                                                      a*f*F*Lr*r*(expl(Power(a,2)*B + S) + expl(Power(a,2)*B + A*Power(r,2) + S) + expl(A*Power(r,2) + T) + 
                                                                  2*expl(Power(a,2)*B + A*Power(r,2) + T))) + 
                                                 L*(Power(a,3)*Power(F,2)*ga*Power(r,2) - Power(a,3)*Power(F,2)*ga*Power(r,2)*expl(Power(a,2)*B) - 
                                                    4*Power(a,2)*f*(-(Fr*r) + F*(-1 + Power(a,2)*Br*r + r*Sr - expl(A*Power(r,2))))*expl(Power(a,2)*B + S) + 
                                                    4*Power(a,2)*F*fr*r*expl(Power(a,2)*B + A*Power(r,2) + S) - 32*expl(Power(a,2)*B + A*Power(r,2) + S + T) + 
                                                    16*Power(a,2)*B*expl(Power(a,2)*B + A*Power(r,2) + S + T) + 
                                                    8*Power(a,3)*Ba*expl(Power(a,2)*B + A*Power(r,2) + S + T) - 
                                                    8*a*Aa*Power(r,2)*expl(Power(a,2)*B + A*Power(r,2) + S + T) - 8*a*Sa*expl(Power(a,2)*B + A*Power(r,2) + S + T) - 
                                                    16*a*Ta*expl(Power(a,2)*B + A*Power(r,2) + S + T) + 32*expl(2*Power(a,2)*B + A*Power(r,2) + S + T))) - 
                                128*a*ga*L*(-1 + expl(Power(a,2)*B))*expl(Power(a,2)*B + A*Power(r,2) + 2*(S + T))))/
            (8.*a*g*L*(-(Power(a,2)*Power(F,2)*g*Power(r,2)) + 16*expl(Power(a,2)*B + A*Power(r,2) + S + T))) ;
            
            
            phival = (-2*a*f*F*g*ha*hr*r + 4*Power(f,2)*Power(hr,2)*expl(Power(a,2)*B + S) + 4*g*Power(ha,2)*expl(A*Power(r,2) + T))/
            (-(Power(a,2)*Power(F,2)*g*L*Power(r,2)) + 16*L*expl(Power(a,2)*B + A*Power(r,2) + S + T)) ;		
            
            /* Nothing should depend on theta once the contractions are performed. */
            
            theta = 0.123456789 ;
            
            
#include "Curv.c"
            
            phi[ii][jj] = phival ;
            
            Riccival = 0. ;
            for(i1=0;i1<5;i1++) {
                for(s1=0;s1<5;s1++) {
                    Riccival += R2[i1][s1]*gu[i1][s1] ;
                }
            }
            
            R2sqrval = 0. ;
            for(i1=0;i1<5;i1++) {
                for(i2=0;i2<5;i2++) {
                    if(R2[i1][i2]!=0.) {
                        for(s1=0;s1<5;s1++) {
                            for(s2=0;s2<5;s2++) {
                                R2sqrval += R2[i1][i2]*R2[s1][s2]*gu[i1][s1]*gu[i2][s2] ;
                            }
                        }
                    }
                }
            }
            
            R4sqrval = 0. ;
            
                
                for(i1=0;i1<5;i1++) {
                    for(i2=0;i2<5;i2++) {
                        for(i3=0;i3<5;i3++) {
                            for(i4=0;i4<5;i4++) {
                                if(C4[i1][i2][i3][i4]!=0.) {
                                    for(s1=0;s1<5;s1++) {
                                        for(s2=0;s2<5;s2++) {
                                            for(s3=0;s3<5;s3++) {
                                                for(s4=0;s4<5;s4++) {
                                                    R4sqrval += C4[i1][i2][i3][i4]*C4[s1][s2][s3][s4]*gu[i1][s1]*gu[i2][s2]*gu[i3][s3]*gu[i4][s4] ;
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
                
           
            
            
            Ricci[ii][jj] = Riccival ;
            
            R2sqr[ii][jj] = R2sqrval ;
            
            R4sqr[ii][jj] = R4sqrval ;
            
            if( phival > phimax ) phimax = phival ;
			
			
		}
	}
    
    
	return( phimax ) ;
	
}






void getG(void)
{
	int ii, jj, var ;
	long pos ;
	
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			for(var=0;var<5;var++) {
				
				if( (pos = mask[var][ii][jj]) == -1 ) continue ;
				
				vecG[pos] = gmet[var][ii][jj] ;
				
			}
			
		}
	}	
	
}


void updateG(REAL mult, REAL *vecptr)
{
	long pos ;
	
	for(pos=0;pos<vecsize;pos++) {
		
		vecG[pos] += mult*vecptr[pos] ;
		
	}	
	
}



void storeG(void)
{
	long pos ;
	int var,ii,jj ;
	
	for(pos=0;pos<vecsize;pos++) {
		
		var   = vecpos[pos][0] ; 
		ii    = vecpos[pos][1] ; 
		jj    = vecpos[pos][2] ; 
		
		gmet[var][ii][jj] = vecG[pos] ;
					
	}
	
}






void init(void)
{
	int ii, jj, var ;
	    
    for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			for(var=0;var<5;var++) { 
				gmet[var][ii][jj] = 0. ; 
                RF[var][ii][jj]   = 0. ;
			}
			
		}
	}
    
}



void initpatches(void)
{
	int ii, jj, var, chk ;
	int intbuf ;
	REAL buf ;
	
	FILE *file ;
	
	
	vecsize = 0 ;
	
	
	file = fopen("patch.dat","r") ;
	
	if(file==NULL) { printf("\n\nError opening patch\n\n") ; exit(1) ; }
	
	fread(&NR, sizeof(int), 1, file) ; fread(&NA, sizeof(int), 1, file) ; 
	
	if( NR>NMAX || NA>NMAX ) { printf("\n\nPatch too large!\n\n") ; exit(1) ; }
	
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			chk = 0 ;
            
            /*
             
             Do not load in coords... but advance the file anyway.
             
             */
            
			chk += fread(&buf, sizeof(double), 1, file) ; 
			chk += fread(&buf, sizeof(double), 1, file) ; 
            
			for(var=0;var<5;var++) {
				chk += fread(&intbuf, sizeof(int), 1, file) ; mask[var][ii][jj] = intbuf ;
				if(intbuf!=-1) {
					if(intbuf>=NMAX2) {
                        
                        printf(" %d \n\n",intbuf) ;
                        
						printf("Vector size too large! %d\n\n",intbuf) ; exit(1) ;
					}
					vecpos[intbuf][0] = var ;
					vecpos[intbuf][1] = ii ;
					vecpos[intbuf][2] = jj ;
					if(intbuf>vecsize) vecsize = intbuf ;
				}
			}
			if(chk!=7) { printf("\n\nError in patch file - %d, %d, %d\n\n",ii,jj,chk) ; exit(1) ; }
			
		}
	}
	
	fclose(file) ;
	
	vecsize += 1 ;
	
	printf("Vector size = %ld\n",vecsize) ;
		
	
}




void save(int writelich)
{
	int ii,jj,var ;
    long posii, posjj, posLich ;
	double buf, val ;
    
    FILE *fileptrLichRow, *fileptrLichCol, *fileptrLichVal ;
		
    
	for(var=0;var<5;var++) {
		for(ii=0;ii<NR;ii++) {
			for(jj=0;jj<NA;jj++) {
                
				buf = (double) gmet[var][ii][jj] ;
				fwrite(&buf, sizeof(double), 1, fileptr) ;	
				
			}
		}
	}
	
	fflush(fileptr) ; 	
	
	for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			buf = (double) phi[ii][jj] ;
			fwrite(&buf, sizeof(double), 1, fileptrphi) ;	
			
		}
	}
	
	fflush(fileptrphi) ; 	
    
    for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			buf = (double) Ricci[ii][jj] ;
			fwrite(&buf, sizeof(double), 1, fileptrRicci) ;	
			
		}
	}
	
	fflush(fileptrRicci) ; 
    
    for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			buf = (double) R2sqr[ii][jj] ;
			fwrite(&buf, sizeof(double), 1, fileptrR2sqr) ;	
			
		}
	}
	
	fflush(fileptrR2sqr) ; 
    
    for(ii=0;ii<NR;ii++) {
		for(jj=0;jj<NA;jj++) {
			
			buf = (double) R4sqr[ii][jj] ;
			fwrite(&buf, sizeof(double), 1, fileptrR4sqr) ;	
			
		}
	}
	
	fflush(fileptrR4sqr) ; 
    
    
    for(var=0;var<5;var++) {
		for(ii=0;ii<NR;ii++) {
			for(jj=0;jj<NA;jj++) {
                
				buf = (double) dvr[var][ii][jj] ;
				fwrite(&buf, sizeof(double), 1, fileptrDerivs) ;	
				
			}
		}
	} 

    for(var=0;var<5;var++) {
		for(ii=0;ii<NR;ii++) {
			for(jj=0;jj<NA;jj++) {
                
				buf = (double) dva[var][ii][jj] ;
				fwrite(&buf, sizeof(double), 1, fileptrDerivs) ;	
				
			}
		}
	} 

    
    for(var=0;var<5;var++) {
		for(ii=0;ii<NR;ii++) {
			for(jj=0;jj<NA;jj++) {
                
				buf = (double) dvrr[var][ii][jj] ;
				fwrite(&buf, sizeof(double), 1, fileptrDerivs) ;	
				
			}
		}
	} 
    
    for(var=0;var<5;var++) {
		for(ii=0;ii<NR;ii++) {
			for(jj=0;jj<NA;jj++) {
                
				buf = (double) dvaa[var][ii][jj] ;
				fwrite(&buf, sizeof(double), 1, fileptrDerivs) ;	
				
			}
		}
	} 
    
    for(var=0;var<5;var++) {
		for(ii=0;ii<NR;ii++) {
			for(jj=0;jj<NA;jj++) {
                
				buf = (double) dvra[var][ii][jj] ;
				fwrite(&buf, sizeof(double), 1, fileptrDerivs) ;	
				
			}
		}
	} 
    
    
	fflush(fileptrDerivs) ; 
    
    
    
    if(writelich==1) {
        
        fileptrLichRow = fopen("LichRow.out","w") ;
        fileptrLichCol = fopen("LichCol.out","w") ;
        fileptrLichVal = fopen("LichVal.out","w") ;
        
        for(posLich=0;posLich<sizeLich;posLich++) {
            
            posii = LichRow[posLich] ;
            posjj = LichCol[posLich] ;
            val   = Lich[posLich] ;
            
            fwrite(&posii, sizeof(long), 1, fileptrLichRow) ;	
            fwrite(&posjj, sizeof(long), 1, fileptrLichCol) ;	
            fwrite(&val, sizeof(double), 1, fileptrLichVal) ;	
            
        }
        
        
        fclose(fileptrLichRow) ; fclose(fileptrLichCol) ; fclose(fileptrLichVal) ;
        
    }
    
}





REAL norm(REAL *vecin)
{
	int ii; 
	REAL tmpval, val ;
	
	for(val=0.0,ii=0;ii<vecsize;ii++) {
		tmpval = fabs(vecin[ii]) ;
		if(val<tmpval) val = tmpval ;
	}
	
	return(val) ;
}


void unpackLich(void)
{
    long posii, posjj, posLich ;
    REAL val ;
     
    for(posii=0;posii<vecsize;posii++) {
        for(posjj=0;posjj<vecsize+1;posjj++) {
            
            
#if LAPACK == 0
            
            Lmat[ posii*(vecsize+1) + posjj ] = 0. ;
                        
#else
            
            Lmat[ posjj*vecsize + posii ] = 0. ;
            
#endif
                
        }        
    }
    
	for(posLich=0;posLich<sizeLich;posLich++) {
		
		posii = LichRow[posLich] ;
		posjj = LichCol[posLich] ;
		val   = Lich[posLich] ;
		
#if LAPACK == 0
        
        Lmat[ posii*(1+vecsize) + posjj ] = val ;
        
#else
        
        Lmat[ posjj*vecsize + posii ] = val ;
        
#endif
        
    }

}



/* LAPACK solver */

void linsolveLAPACK(REAL *x, REAL *b)
{
    int N, nrhs, ldLmat, ldx, info ;
    int ipiv[NMAX2];
    
    long ii ;
    
    N    = vecsize ;
    nrhs = 1 ;
    ldLmat  = vecsize ;
    ldx  = vecsize ;
    
    for(ii=0;ii<vecsize;ii++) x[ii] = b[ii] ;
    
    dgesv_(&N, &nrhs, Lmat, &ldLmat, ipiv, x, &ldx, &info);
    
    if(info != 0) {
        
        printf("\n\nFailure to solve linear system - %d\n\n", info);
        
        exit(1) ;
        
    }
    
    
}


/* Full/Partial pivoting */

void linsolveDirect(REAL *x, REAL *b)
{
    long ii, jj, kk, iimax, jjmax ;
    
    int tmp, index[NMAX2] ;
    
    REAL maxval, val, y[NMAX2] ;
    
    REAL tmpval, tmp2, tmp3 ;
    long pos1, pos2 ;
    
    
    for(ii=0;ii<vecsize;ii++) { index[ii] = ii ; }
    
    for(ii=0;ii<vecsize;ii++) { Lmat[ ii*(vecsize+1) + vecsize ] = b[ii] ; } 
    

  
    
    printf("\n<") ;
    for(kk=0;kk<vecsize;kk++) if( kk%(vecsize/100)==0 ) printf(" ") ;
    printf(">\n") ;
    
    
    
    printf("\n<") ;
	for(kk=0;kk<vecsize;kk++) {
        if( kk%(vecsize/100)==0) printf(".") ; fflush(stdout) ;
        
    
        
// Partial pivot on rows
        
        maxval = 0 ; iimax = -1 ;
        
        for(ii=kk;ii<vecsize;ii++) {
            
            if( fabs( Lmat[ii*(vecsize+1)+kk] ) > maxval ) {
                
                maxval = fabs(Lmat[ii*(vecsize+1)+kk]) ;
                
                iimax = ii ;
                
            }
            
        }
        
        if( iimax == -1 ) {
            
            printf("\n\nSingular Lich matrix!\n\n") ;
            exit(1) ;
            
        }
        
// Swap rows

        for(jj=kk;jj<vecsize+1;jj++) {
            
            val = Lmat[iimax*(vecsize+1)+jj] ;
            
            Lmat[iimax*(vecsize+1)+jj] = Lmat[kk*(vecsize+1)+jj] ;
            
            Lmat[kk*(vecsize+1)+jj] = val ;
            
        }
        

        
        
        pos2 = kk*(vecsize+1)+kk+1 ;
 
        tmp3 = -1./Lmat[kk*(vecsize+1)+kk] ;
        
        for(ii=kk+1;ii<vecsize;ii++) {
            
            tmpval = Lmat[ii*(vecsize+1)+kk]*tmp3 ;
            
            pos1 = ii*(vecsize+1)+kk+1 ;
            
            if(tmpval==0.) continue ;
            
            for(jj=0;jj<vecsize-kk;jj++) {
                
                tmp2 = Lmat[pos2+jj] ;
                
             // Slower!!  
                //if(tmp2==0.) continue ;
                
                Lmat[pos1+jj] += tmp2*tmpval ;
                
            }
            
            Lmat[ii*(vecsize+1)+kk] = 0. ;
            
        }
        
        
    }
    
    
    printf(">\n") ;
    
    
// Back substitution
    
    
    for(kk=vecsize-1;kk>=0;kk--) {
    
        y[kk] = Lmat[ kk*(vecsize+1)+vecsize ]   ;
        
        for(ii=kk+1;ii<vecsize;ii++) {
        
            y[kk] += - Lmat[kk*(vecsize+1)+ii]*y[ii] ; 
            
        }
        
        y[kk] /= Lmat[kk*(vecsize+1)+kk] ;
        
        x[index[kk]] = y[kk] ;

    }
    
}


// Wrapper

void linsolve(REAL *x, REAL *b)
{
    
#if LAPACK == 0
    
    linsolveDirect(x, b) ; 
    
#else
    
    linsolveLAPACK(x, b) ; 
    
#endif   
    
}



void linsolvechk(REAL *x, REAL *b)
{
    long ii, jj ;
    
    REAL val, y[NMAX2] ;

    
// Check accuracy
    
    val = 0 ;
    
    for(ii=0;ii<vecsize;ii++) {
        
        y[ii] = b[ii] ;
        
        for(jj=0;jj<vecsize;jj++) {
            
            
            #if LAPACK == 0
         
            y[ii] += - Lmat[ii*(vecsize+1)+jj]*x[jj] ; 
            
            #else
            
            y[ii] += - Lmat[jj*vecsize+ii]*x[jj] ; 
            
            #endif
            
        }
        
        if( fabs(y[ii]) > val ) val = fabs(y[ii]) ;
        
    }
    
    
#if REALTYPE == 2
    
    printf("\n   Lin solve check ->  %10.3Le\n",val) ;
    
    
#else   
    
    printf("\n   Lin solve check ->  %10.3e\n",val) ;
    
   
#endif 

	
}




void setupCheb(void)
{
    int Neven ;
    
    int ii, jj, kk ;
    REAL x[NMAX], c1, c2 ;
    
    REAL Dmat[NMAX][NMAX], D2mat[NMAX][NMAX] ;
    
    
    /* Std - for r  */
    
    for(ii=0;ii<NR;ii++) {
        
        x[ii] = (1 + cosl(PI*(((REAL) (NR-ii))-1)/(((REAL) NR)-1)) )/2 ;
        
    }
    
    for(ii=1;ii<NR-1;ii++) {
		for(jj=1;jj<NR-1;jj++) {
            
            if(ii==jj) 
                Dmat[ii][jj] = - (x[jj] - 0.5)/(2*(1-x[jj])*x[jj]) ;
            else
                Dmat[ii][jj] = ((REAL) pow(-1,ii+jj))/( x[ii] - x[jj] ) ;
            
        }
    }
    
    Dmat[0][0] = - ( 2*((REAL) pow(NR-1,2)) + 1 )/3 ;
    
    for(jj=1;jj<NR-1;jj++) {
        
        Dmat[0][jj] = -2*((REAL) pow(-1.,jj))/x[jj] ;
        
    }
    
    Dmat[0][NR-1] = - pow(-1,NR-1) ;
    
    Dmat[NR-1][NR-1] = + ( 2*((REAL) pow(NR-1,2)) + 1 )/3 ;
    
    for(jj=1;jj<NR-1;jj++) {
        
        Dmat[NR-1][jj] = 2*((REAL) pow(-1,NR+jj-1))/(1 - x[jj]) ;
        
    }
    
    Dmat[NR-1][0] = + pow(-1,NR-1) ;
    
    for(ii=1;ii<NR-1;ii++) {
        
        Dmat[ii][0] = +((REAL) pow(-1,ii))/x[ii]/2 ;
        
    }
    
    for(ii=1;ii<NR-1;ii++) {
        
        Dmat[ii][NR-1] = - ((REAL) pow(-1,NR+ii-1))/(1 - x[ii])/2 ;
        
    }
    
    for(ii=0;ii<NR;ii++) {
        for(jj=0;jj<NR;jj++) {
            
            D2mat[ii][jj] = 0. ;
            
            for(kk=0;kk<NR;kk++) {
                D2mat[ii][jj] += Dmat[ii][kk]*Dmat[kk][jj] ;
            }
            
        }
    }
    
    
    for(ii=0;ii<NR;ii++) {
        
		for(jj=0;jj<NR;jj++) {
            
#if EvenR != 1            
            
            patch[ii][jj][0] = x[ii] ;
            
            DRmat[ii][jj] = Dmat[ii][jj] ; 
            
            DR2mat[ii][jj] = D2mat[ii][jj] ; 
            
#endif              
            
            
#if EvenA != 1
            
            patch[ii][jj][1] = x[jj] ;
            
            DAmat[ii][jj] = Dmat[ii][jj] ; 
            
            DA2mat[ii][jj] = D2mat[ii][jj] ;
            
#endif                
            
            
        }
    }
    
    
    
    
    /* Even for a */
    
    
    Neven = 2*NR - 1 ;
    
    for(ii=0;ii<=Neven;ii++) {
        
        x[ii] = cosl(PI*(((REAL) (Neven-ii)))/Neven) ;
        
    }
    
    for(ii=0;ii<=Neven;ii++) {
		for(jj=0;jj<=Neven;jj++) {
            
            if(ii==0 || ii==Neven) c1=2 ; else c1=1 ;
            
            if(jj==0 || jj==Neven) c2=2 ; else c2=1 ;
            
            if(ii!=jj) Dmat[ii][jj] = ((REAL) c1)/((REAL) c2)*((REAL) pow(-1,ii+jj))/( x[ii] - x[jj] ) ;
            
        }
    }
    
    for(ii=1;ii<=Neven-1;ii++) Dmat[ii][ii] = - x[ii]/(2*(1-x[ii]*x[ii])) ; 
    
    Dmat[0][0] = - ( 2*((REAL) pow(Neven,2)) + 1 )/6 ;
    
    Dmat[Neven][Neven] = + ( 2*((REAL) pow(Neven,2)) + 1 )/6 ;
    
    for(ii=0;ii<=Neven;ii++) {
        for(jj=0;jj<=Neven;jj++) {
            
            D2mat[ii][jj] = 0. ;
            
            for(kk=0;kk<=Neven;kk++) {
                D2mat[ii][jj] += Dmat[ii][kk]*Dmat[kk][jj] ;
            }
            
        }
    }
    
    
    
    
    for(ii=0;ii<NR;ii++) {
        
		for(jj=0;jj<NR;jj++) {
            
#if EvenR == 1
            
            patch[ii][jj][0] = x[ii+NR] ;
            
            DRmat[ii][jj] = Dmat[ii+NR][jj+NR] + Dmat[ii+NR][NR-1-jj] ; 
            
            DR2mat[ii][jj] = D2mat[ii+NR][jj+NR] + D2mat[ii+NR][NR-1-jj] ;
            
#endif             
            
            
#if EvenA == 1
            
            patch[ii][jj][1] = x[jj+NR] ;
            
            DAmat[ii][jj] = Dmat[ii+NR][jj+NR] + Dmat[ii+NR][NR-1-jj] ; 
            
            DA2mat[ii][jj] = D2mat[ii+NR][jj+NR] + D2mat[ii+NR][NR-1-jj] ;
            
#endif            
            
        }
    }
    
    
    
    
}