Programa en C

/*
   UNIVERSIDAD AUTONOMA DE PUEBLA
       COLEGIO DE COMPUTACION
        ANALISIS NUMERICO II
*/

#include <stdio.h>
#include <math.h>
#include <ctype.h>
#include "plotp.h"
#include "arcom.c"

#define NF  	29
#define NC  	29

complex inc, z0;
double XMI, XMA, IMI, IMA, DMIx, DMAx, DMIi, DMAi, N, NDIVX, NDIVI;

 void ejesri()
 { pltig( XMI, IMI, 1, pltca);
   pltig( XMA, IMA, 2, pltca);
   pltig( XMI, 0.0, 3, pltca);
   pltig( XMA, 0.0, 8, pltca);
   pltig( 0.0, IMI, 3, pltca);
   pltig( 0.0, IMA, 9, pltca);
   return;  }

complex func(z) complex z; {                            /*  Calculo de la funci\'on   */
complex fz={0.0,0.0},caux,ca={1.0,0.0},com;
com.x=DMIx;  com.i=DMIi;
  for(com.x=DMIx; com.x<=DMAx; com.x+=inc.x)
    for(com.i=DMIi; com.i<=DMAi; com.i+=inc.i) {
      caux=dfnc(dfnc(z,com),z0);
      fz=adnc(fz,divnc(ca,mulnc(caux,caux))); };
/*  return(divnc(ca,fz)); */                             /* 1/f(z)                     */
  return(fz);                                            /* f(z)                       */
com.x=1; com.i=1;
/*return(mulnc(fz, mulnc(dfnc(z, com),dfnc(z, com))));*/ /*factorizando el polo        */
/*return(mulnc(fz, dfnc(z, com)));                    */ /*factorizando un orden del polo*/
/*com.x=0.5; com.i=0.5;
  return(divnc(fz, mulnc(dfnc(z,com), dfnc(z,com))));*/ /* factoriznado un cero      */
/*return(divnc(fz, dfnc(z,com)));	             */ /*factorizando un orden del cero*/
  }


double acota(double n)
{
  if (n>N) return (N);
  else if( n<-N ) return (-N);
  return (n);
}

void main()
{ int i,j,decis;
  complex aux;
  double mfdp[NC][NF], mcf[NC][NF];
  double dx,di,zmi,zma,x,ic,dz;
  double micx, macx, mici, maci, smicx, smacx, smici, smaci, sdct;

printf("Proporcionar XMI, XMA, IMI, IMA, z0 :   ");
scanf("%lf %lf %lf %lf %lf %lf", &XMI, &XMA, &IMI, &IMA, &z0.x, &z0.i);
printf("Proporcionar Cota, NDIVX, NDIVI, DMIx, DMAx, DMIi, DMAi: ");
scanf("%lf %lf %lf %lf %lf %lf %lf", &N, &NDIVX, &NDIVI, &DMIx, &DMAx, &DMIi, &DMAi);

printf("Rango de primer corte :  ");
scanf("%lf %lf %lf %lf", &micx, &macx, &mici, &maci);

printf("Rango de segundo corte :  ");
scanf("%lf %lf %lf %lf", &smicx, &smacx, &smici, &smaci);

inc.x=(DMAx-DMIx)/(NDIVX); inc.i=(DMAi-DMIi)/(NDIVI);
  dx=fabs(XMA-XMI)/(double)(NF-1);
  di=fabs(IMA-IMI)/(double)(NC-1);
  zmi=acota(modnc(func(XMI,IMI)));
  zmi=(zmi<0.0) ? 0.0 : zmi;
  zma = zmi;

  sdct=N/4;

  ic = IMI;
  for(i=0; i<NC; i+=1)
  {
    x = XMI;
    for(j=0; j<NF; j+=1) {
    if((micx<=x && x<=macx) && (mici<=ic && ic<=maci)) {   /*Para primer corte        */
      mfdp[i][j]=-1;
      x+=dx;
      continue;  };
    aux=arg_mod(func(x,ic));
    mcf[i][j]=aux.i;
    aux.x=x;   aux.i=ic;
    aux.x=acota(modnc(func(aux)));
    if((smicx<=x&&x<smacx) && (smici<=ic&&ic<smaci))       /*Para segundo corte       */
      mfdp[i][j] = (aux.x>sdct)?sdct:aux.x;
    else
      mfdp[i][j] = aux.x;

    if (mfdp[i][j] < zmi) zmi=mfdp[i][j];
    if (mfdp[i][j] > zma) zma=mfdp[i][j];
    x+=dx;   }
    ic+=di;
   }

if(toupper(decis)=='S')  {
for(i=0;i<NC;i+=NC-1)                                      /*Para los bordes          */
  for(j=0;j<NF;mfdp[i][j++]=zmi);
for(j=0;j<NF;j+=NF-1)
  for(i=0;i<NC;mfdp[i++][j]=zmi);  }

for(i=0;i<NC;i++)                                          /*Para primer corte        */
  for(j=0;j<NF;j++) {
    if(mfdp[i][j]!=-1)
      continue;
    mfdp[i][j]=zmi;  };

   dz = fabs(zma-zmi)/8.0;
   x = zmi - dz;
   i = zma + dz;

decis = 'a';
while(decis != 'X') {
printf("\n\t-S- Superficie\n\t-F- cont. Fase\n\t-2- 2contorno\n\t-C- Contorno");
decis=toupper(getche());
   pltdk(3);
   plt00();
   plotf(-1);
   pltfr(3);
   pltso(3.0, 3.0, 0.0, 0.0);
   switch (decis) {
     case 'C':
          ejesri();
          pltcolor(14);                                /*  YELLOW                    */
          pltkp(zmi,&mfdp[0][0],(zma-zmi)/150,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
          pltcolor(8);                                 /*  DARKGRAY                  */
          pltkp((zma-zmi)/150,&mfdp[0][0],(zma-zmi)/2,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
          pltcolor(9);                                 /* BROWN                      */
          pltkp((zma-zmi)/2,&mfdp[0][0],2*(zma-zmi)/3,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
          pltcolor(1 | 0);                             /* BLUE | BLACK               */
          pltkp(2*(zma-zmi)/3,&mfdp[0][0],zma,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
	       
            /* Para graficar la red de periodos  */
          pltcolor(4);                                 /* RED                        */
          DMAx=(DMAx>XMA)?DMAx:XMA;
          DMAi=(DMAi>IMA)?DMAi:IMA;
          DMIx=(DMIx<XMI)?DMIx:XMI;
          DMIi=(DMIi<IMI)?DMIi:IMI;
          for(dx=DMIx+(inc.x/2); dx<=DMAx; dx+=inc.x)  {
            pltig(dx,DMIi, 3, pltca);
            pltig(dx,DMAi, 4, pltca);         }
          for(di=DMIi+(inc.x/2); di<=DMAi; di+=inc.i)  {
            pltig(DMIx, di, 3, pltca);
	    pltig(DMAi, di, 4, pltca);        }
            pltcolor(11);
	    ejesri();
            break;
     case 'F': 
          ejesri();
          pltkc(zmi,&mcf[0][0],zma,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
          break;
     case '2': 
          ejesri();
	  pltcolor(14);                                /* YELLOW                     */
          pltkp(zmi,&mfdp[0][0],(zma-zmi)/150,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
          pltcolor(8);                                 /* DARKGRAY                   */
          pltkp((zma-zmi)/150,&mfdp[0][0],(zma-zmi)/2,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
	  pltcolor(9);                                 /* BROWN                      */
          pltkp((zma-zmi)/2,&mfdp[0][0],2*(zma-zmi)/3,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
	  pltcolor(1 | 0);                             /* BLUE | BLACK               */
          pltkp(2*(zma-zmi)/3,&mfdp[0][0],zma,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
          pltcolor(2);                                 /* GREEN                      */
	  pltkc(zmi,&mcf[0][0],zma,25,(NF/64)+1,NF,(NC/64)+1,NC,pltca);
	  pltcolor(11);
	  ejesri();
	  break;
     case 'S': 
          for(i=0;i<NC;i+=NC-1)                        /* Para los bordes            */
	     for(j=0;j<NF;mfdp[i][j++]=zmi);
	     for(j=0;j<NF;j+=NF-1)
		for(i=0;i<NC;mfdp[i++][j]=zmi);
          pltrv(x,&mfdp[0][0], i, NF, NC, -40.0, pltca);
          break;
     case 'X': 
          exit(0);
     default : 
          printf("No detecto");
          break;  };

   getchar();
   pltcolor(2);  
   plotr();
   pltla(" fndp ",6);
   pltho();
   closeplt(); }
   }