home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / trkman_map.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-29  |  4.3 KB  |  162 lines
  1. /*
  2. ----------------------------------------------------------------------
  3.      TRANLATION OF TRKMANQ.F INTO C 
  4.      Compute only one set of manifolds and create a copy by iteration
  5. ----------------------------------------------------------------------
  6.      compute manifolds of ANY given set of hyperbolic fixed points
  7.      with eigenvalues and eigenvectors
  8.      m[us]f : finer divisions of manifolds
  9.      iskip : number of skips between each data
  10. ----------------------------------------------------------------------
  11. NOTE: In case of a negative eigenvalue use rnew = r * 2 
  12. ----------------------------------------------------------------------
  13. */
  14. #include "../include/x11r2_kaos_def.h"
  15.  
  16. #define DELFRAC 0.1
  17. void trkman_map(sevec,seval,xe,r,ms,mu,msf,muf,iskip,delman,n)
  18. int r,ms,mu,msf,muf,iskip,n;
  19. double **sevec,**seval,*xe,delman;
  20. {
  21.     int i,ic,jc,mc,kc,mf,man_index,icnt,color_index,color,ndid,rnew,negative_ev_on;
  22.     double fabs(),delx,xerr,factor0,factor,**x,**xp,*x1,*x2,exp(),log(),*dvector(),**dmatrix();
  23.     extern int stop,dot_size,fp_display_option,fi_maxsq;
  24.     extern double fi_eps,fi_epsf,fi_epsm;
  25.     extern char string[];
  26.  
  27.     x1 = dvector(0,n-1);
  28.     x2 = dvector(0,n-1);
  29.     mf = (muf > msf ? muf : msf);
  30.     x = dmatrix(0,mf-1,0,n-1);
  31.     xp = dmatrix(0,mf-1,0,n-1);
  32.  
  33.     for(ic=0;ic<n;ic++){    
  34.         /* if complex eigenvalue exit*/
  35.         if(seval[ic][1]!=0)
  36.             man_index = -1;
  37.         else {
  38.             if(fabs(seval[ic][0])>1){
  39.                 man_index = 0;
  40.                 color = Blue;
  41.             }
  42.             else {
  43.                 man_index = 1;
  44.                 color = Red;
  45.             }
  46.             if(seval[ic][0]<0){
  47.                 negative_ev_on=1;
  48.                 rnew = r * 2;
  49.             }
  50.             else {
  51.                 negative_ev_on=0;
  52.                 rnew = r;
  53.             }    
  54.         }
  55.         if(man_index==0) {
  56.             if(negative_ev_on){
  57.                 factor0=exp(log(seval[ic][0])*2./muf);
  58.             }
  59.             else{
  60.                 factor0=exp(log(seval[ic][0])/muf);
  61.             }
  62.         }
  63.         else if(man_index==1){
  64.             if(negative_ev_on){
  65.                 factor0=exp(log(seval[ic][0])*2./msf);
  66.             }
  67.             else{
  68.                 factor0=exp(log(seval[ic][0])/msf);
  69.             }
  70.         }
  71.  
  72.         sprintf(string,"Computing a manifold for ev[%d]...",ic);
  73.         printf("%s\n",string);
  74.         if(man_index==0){
  75.             printf("(Ev=%g %g)\n",seval[ic][0],seval[ic][1]);
  76.             printf("(Evec %d =%g %g)\n",ic,sevec[ic][0],sevec[ic][1]);
  77.             for(jc=0;jc<2;jc++){
  78.                 factor = 1.;
  79.                 for(mc=0;mc<muf;mc++){
  80.                     if(jc==0){
  81.                         printf("Unstable manfold for ev[%d]...\n",ic);
  82.                         for(i=0;i<n;i++) x[mc][i] = xe[i]+delman*sevec[ic][i]*factor;
  83.                     }
  84.                     else {
  85.                         printf("Unstable manfold for ev[%d]...\n",ic);
  86.                         for(i=0;i<n;i++) x[mc][i] = xe[i]-delman*sevec[ic][i]*factor;
  87.                     }
  88.                     factor *= factor0;
  89.                 }
  90.                 icnt=0;
  91.                 for(kc=0;kc<mu;kc++){
  92.                     for(mc=0;mc<muf;mc++){
  93.                         fmap_user(x[mc],rnew,n);
  94.                         for(i=0;i<n;i++)x1[i]=x[mc][i];
  95.                         if((icnt % iskip)==0){
  96.                             (void) draw_record_pwf(x1,color,-1,dot_size,0,1);
  97.                               if(fp_display_option==1)
  98.                                   (void) draw_record_other_pwf(x1,rnew,color,-1,dot_size,0,1);
  99.                         }
  100.                         icnt++;
  101.                         if(stop){
  102.                             goto done;
  103.                         }
  104.                     }
  105.                 }
  106.             }
  107.         }
  108.         else if(man_index==1) {
  109.             printf("(Ev=%g %g)\n",seval[ic][0],seval[ic][1]);
  110.             printf("(Evec %d =%g %g)\n",ic,sevec[ic][0],sevec[ic][1]);
  111.             for(jc=0;jc<2;jc++){
  112.                 icnt=0;
  113.                 factor=1;
  114.                 for(mc=0;mc<msf;mc++){
  115.                     for(i=0;i<n;i++) xp[mc][i]=xe[i];
  116.                     if(jc==0){
  117.                         printf("Stable manfold for ev[%d]...\n",ic);
  118.                         for(i=0;i<n;i++) x[mc][i] = xp[mc][i]+delman*sevec[ic][i]/factor;
  119.                     }
  120.                     else {
  121.                         printf("Stable manfold for ev[%d]...\n",ic);
  122.                         for(i=0;i<n;i++) x[mc][i] = xp[mc][i]-delman*sevec[ic][i]/factor;
  123.                     }
  124.                     factor *= factor0;
  125.                 }
  126.                 for(kc=0;kc<ms;kc++){
  127.                     for(mc=0;mc<msf;mc++){
  128.                         /* install adjustable delc congtroller */
  129.                         /* currently it is divided by the same eigenvalue */
  130.                         for(i=0;i<n;i++){
  131.                             delx=(x[mc][i]-xp[mc][i])/seval[ic][0];
  132.                             x2[i]=x[mc][i]+delx;
  133.                             x1[i]=x1[i]- DELFRAC*delx;
  134.                         }
  135.                         (void) fmapi_user(x1,x2,x[mc],fi_maxsq,n,rnew,fi_eps,fi_epsf,&xerr,&ndid);
  136.                         if(stop){
  137.                             goto done;
  138.                         }
  139.                         for(i=0;i<n;i++){
  140.                             xp[mc][i]=x[mc][i];
  141.                             x[mc][i]=x1[i];
  142.                         }
  143.  
  144.                         if((icnt%iskip)==0){
  145.                             (void) draw_record_pwf(x1,color,-1,dot_size,0,1);
  146.                               if(fp_display_option=1)
  147.                                   (void) draw_record_other_pwf(x1,rnew,color,-1,dot_size,0,1);
  148.                         }
  149.                         icnt++;
  150.                     }
  151.                 }
  152.             }
  153.         }
  154.     }
  155.  
  156. done:
  157.     (void) free_dvector(x1,0,n-1);
  158.     (void) free_dvector(x2,0,n-1);
  159.     (void) free_dmatrix(x,0,mf-1,0,n-1);
  160.     (void) free_dmatrix(xp,0,mf-1,0,n-1);
  161. }
  162.