/*============================= PNP tester ====================================*/#include "PNP0.C"#include "stdlib.h"#include "stdio.h"static void def_comp_TC(PNP_SMP *smp){L kion, k;for (kion=0; kion<smp->Nions; kion++)	for (k=0; k<smp->Ngrid; k++) ((smp->TCs)[kion])[k] = ((smp->Cs)[kion])[k];}static void do_angle(DF *phi, L n){L k; DF pi_2 = Pi / 2;for (k=0; k<n; k++)   { phi[k] = ((DF)k + 1) / ((DF)n + 2) * pi_2;   }}void main(void){PNP_SMP *smp;L sel, k, kf, kl;L nions, ngrid, npore, indices[5], dimensions[5];DF poreradius, porelength, lengths[5];while (true)  { printf("\n# ");    if (scanf("%d",&sel) == 1)     { switch(sel)        {        case 1: nions = 2; ngrid = 1000; poreradius = 3.0; porelength = 10.0;                smp = build_cyl_geom(nions, ngrid, poreradius, porelength,                                      malloc);                printf("\nmem = %d, %d, %d",smp,(L)smp->memory,smp->nbytes);                smp->comp_TC = def_comp_TC; break;        case 2: (smp->Z)[0] = -1.0; (smp->Z)[1] = 1.0;                (smp->CL)[0] = 0.3; (smp->CL)[1] = 0.3;                (smp->CR)[0] = 0.3; (smp->CR)[1] = 0.3;                smp->VM = -0.1; smp->Celsius = 10.0;                for (k=0; k<ngrid; k++)                  { (smp->F)[k] = -5.0;                    (smp->EPS)[k] = 80.0;                    ((smp->Ds)[0])[k] = 2e-8;                    ((smp->Ds)[1])[k] = 1e-8;                    ((smp->E0s)[0])[k] = 0.0;                    ((smp->E0s)[1])[k] = 0.0;                  }                 break;        case 3: smp->MaxIter = 50; smp->Tolerance = 1e-10;                if (solve_PNP(smp))                   { printf("\nNetCurrent = %.14Le",smp->NetCurrent);                     printf("\nFluxes = %e, %e",smp->FLUX[0],smp->FLUX[1]);                   } else printf("\nNo convergence!");                break;        case 11: nions = 2; npore = 1000; poreradius = 3.0;                 lengths[0] = 100;                 lengths[1] = 20;                 lengths[2] = 10;                 lengths[3] = 20;                 lengths[4] = 100;                 smp = build_cab_geom(nions, npore, poreradius, lengths, do_angle,                       do_angle, malloc, free, indices, dimensions);                  printf("\nmem = %d, %d, %d",smp,(L)smp->memory,smp->nbytes);                 smp->comp_TC = def_comp_TC;  break;                               case 12: (smp->Z)[0] = -1.0; (smp->Z)[1] = 1.0;                 (smp->CL)[0] = 0.3; (smp->CL)[1] = 0.3;                 (smp->CR)[0] = 0.3; (smp->CR)[1] = 0.3;                 smp->VM = -0.1; smp->Celsius = 10.0;                 for (k=0; k<smp->Ngrid; k++)                  { (smp->F)[k] = 0.0;                    (smp->EPS)[k] = 80.0;                    ((smp->Ds)[0])[k] = 2e-6;                    ((smp->Ds)[1])[k] = 1e-6;                    ((smp->E0s)[0])[k] = 0.0;                    ((smp->E0s)[1])[k] = 0.0;                  }                 kf = indices[2]; kl = kf + dimensions[2];                 for (k=kf; k<kl; k++)                  { (smp->F)[k] = -5.0;                    ((smp->Ds)[0])[k] = 2e-8;                    ((smp->Ds)[1])[k] = 1e-8;                  }                 break;        default: printf("...no!"); break;        }     }  }}