home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / examples / demo / demosrc / d_optimize.pro < prev    next >
Encoding:
Text File  |  1997-07-08  |  89.3 KB  |  2,962 lines
  1. ; $Id: d_optimize.pro,v 1.17 1997/04/21 19:52:51 tremblay Exp $
  2. ;
  3. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;
  6. ;+
  7. ;  FILE:
  8. ;       d_optimize.pro
  9. ;
  10. ;  CALLING SEQUENCE: d_optimize
  11. ;
  12. ;  PURPOSE:
  13. ;       Shows the path taken by 2 optimization algorithms :
  14. ;       Powell and DFP.
  15. ; 
  16. ;  MAJOR TOPICS: Data analysis and plotting
  17. ;
  18. ;  CATEGORY:
  19. ;       IDL 5.0
  20. ;
  21. ;  INTERNAL FUNCTIONS and PROCEDURES:
  22. ;       pro opt_init16Colors    - Initialize 16 working colors
  23. ;       pro opt_initPoint       - Initialize the starting points coordinates
  24. ;       pro opt_resetPoint      - Redraw the starting point
  25. ;       pro opt_resetSurface    - Redraw the surface (function)
  26. ;       pro opt_resetResults    - Reset the result labels
  27. ;       pro opt_showStart       - Show the starting point (Label and point)
  28. ;       pro about_Optimization  - Display the inforation file
  29. ;       pro opt_plotPoint       - Plot the points along the path
  30. ;       pro DFP                 - DFP algorithm
  31. ;       pro poweop              - Powell algorithm
  32. ;       pro optimize_Event      - Event handler
  33. ;       pro optimize_Cleanup    - Cleanup
  34. ;       pro d_optimize            - Main procedure
  35. ;
  36. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  37. ;       pro funcop              - function to minimize (DFP and Powell)
  38. ;       pro linmop              - subroutine in Powell
  39. ;       pro breop5              - subroutine in Powell
  40. ;       pro mnbop5              - subroutine in Powell
  41. ;       pro f1dop5              - subroutine in Powell
  42. ;       pro sign                - subroutine in Powell
  43. ;       fun gettips             - Get the tip text structure.
  44. ;       pro widtips             - Create the widget text for tips.
  45. ;       pro widtips             - Size the widget text for tips.
  46. ;       optimize.txt
  47. ;       optimize.tip
  48. ;
  49. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  50. ;             Numerical recipes in C, 2nd ed. 
  51. ;             Press, Vetterling, Teukolsky, and Flannery
  52. ;             Cambridge University Press
  53. ;             ISBN 0-521-43108-5
  54. ;
  55. ;  NAMED STRUCTURES:
  56. ;       none.
  57. ;
  58. ;  COMMON BLOCS:
  59. ;       none.
  60. ;
  61. ;  MODIFICATION HISTORY:  Written by: DAT,RSI, 1996
  62. ;
  63. ;-
  64.  
  65. ; $Id: d_optimize.pro,v 1.17 1997/04/21 19:52:51 tremblay Exp $
  66. ;
  67. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  68. ;       Unauthorized reproduction prohibited.
  69. ;
  70. ;+
  71. ;  FILE:
  72. ;       dfp.pro
  73. ;
  74. ;  CALLING SEQUENCE: dfp
  75. ;
  76. ;  PURPOSE:  Compute the minimum of a given function
  77. ;            using the Dorovan-Fletcher-Powell (DFP) algorithm.
  78. ;            Specific use for optim5.pro
  79. ;
  80. ;  MAJOR TOPICS: Data analysis
  81. ;
  82. ;  CATEGORY:
  83. ;       IDL 4.0
  84. ;
  85. ;  INTERNAL FUNCTIONS and PROCEDURES:
  86. ;       fun grad                - Returns the gradient
  87. ;       pro lnsrch              - Perform the line search
  88. ;       pro dfp                 - Main procedure
  89. ;
  90. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  91. ;       pro funcop              - Returns the funciton value
  92. ;
  93. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  94. ;             Numerical recipes in C, 2nd ed.
  95. ;             Press, Vetterling, Teukolsky, and Flannery
  96. ;             Cambridge University Press
  97. ;             ISBN 0-521-43108-5
  98. ;
  99. ;  NAMED STRUCTURES:
  100. ;       none.
  101. ;
  102. ;  COMMON BLOCS:
  103. ;       none.
  104. ;
  105. ;  MODIFICATION HISTORY:Written by DAT,RSI, June 1996
  106. ;-
  107. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  108. ;       Unauthorized reproduction prohibited.
  109. ;
  110. ;+
  111. ;  FILE:
  112. ;       funcop.pro
  113. ;
  114. ;  CALLING SEQUENCE: funcop
  115. ;                    Called by : dfp.pro
  116. ;                                f1dop5.pro
  117. ;                                lnarch.pro
  118. ;                                poweop.pro
  119. ;
  120. ;  PURPOSE:  Returns the function value.
  121. ;            Specific use for optim5.pro
  122. ;
  123. ;  MAJOR TOPICS: Data analysis
  124. ;
  125. ;  CATEGORY:
  126. ;       IDL 4.0
  127. ;
  128. ;  INTERNAL FUNCTIONS and PROCEDURES:
  129. ;       NONE
  130. ;
  131. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  132. ;       NONE
  133. ;
  134. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  135. ;             Numerical recipes in C, 2nd ed.
  136. ;             Press, Vetterling, Teukolsky, and Flannery
  137. ;             Cambridge University Press
  138. ;             ISBN 0-521-43108-5
  139. ;
  140. ;  NAMED STRUCTURES:
  141. ;       none.
  142. ;
  143. ;  COMMON BLOCS:
  144. ;       none.
  145. ;
  146. ;  MODIFICATION HISTORY:Written by DAT,RSI, June 1996
  147. ;-
  148. ;---------------------------------------------------------------------------
  149. ;
  150. ;  Purpose:  Returns the funciton value at point X
  151. ;
  152. function FUNCOP, $
  153.     X, $           ; IN: data coordinates
  154.     surfaceIndex   ; IN: surface index (1, 2, or 3)
  155.  
  156.     x(0) = double(x(0))
  157.     x(1) = double(x(1))
  158.  
  159.     ;  Banana function.
  160.     ;
  161.     if (surfaceIndex EQ 1 ) then begin
  162.         val = 100.0D0*(x(1)-x(0)*x(0))^2 + (1.0-x(0))^2
  163.         RETURN, DOUBLE(val)
  164.     endif
  165.  
  166.     if (surfaceIndex EQ 2 ) then begin
  167.         RETURN, 3.0 * (1.0-x(0))^2 * EXP(-x(0)^2 - (x(1)+1.0)^2)            $
  168.             -10.0 * (x(0)/5.0 - x(0)^3 - x(1)^5) * EXP(-x(0)^2 - x(1)^2) $
  169.             -EXP(-(x(0)+1.0)^2 - x(1)^2)/3.0
  170.     endif
  171.  
  172.     if (surfaceIndex EQ 3 ) then begin
  173.         RETURN, DOUBLE( (x(0) + 2.0 *x(1) ) * exp(-x(0)^2 - x(1)^2) ) 
  174.     endif
  175.  
  176. end  ; Of FUNCOP
  177.  
  178.  
  179. ;---------------------------------------------------------------------------
  180. ;
  181. ;  Purpose: Returns the gradient value at point X
  182. ;
  183. function GRAD, $
  184.     X, $          ; IN: data point
  185.     surfaceindex  ; IN: surface index (1, 2, or 3).
  186.  
  187.     x(0) = double(x(0))
  188.     x(1) = double(x(1))
  189.  
  190.     ;  Function # 1  : banana function.
  191.     ;
  192.     if (surfaceindex eq 1) then begin
  193.          dzdx0= -400.0*x(0)*(x(1)-x(0)*x(0)) - 2*(1-x(0))
  194.          dzdx1= 200.0*(x(1)-x(0)*x(0))
  195.          RETURN, [double(dzdx0),double(dzdx1)]
  196.     endif
  197.  
  198.     ;  Function # 2 .
  199.     ;
  200.     if (surfaceindex eq 2) then begin
  201.         exp0 = exp(-x(0)^2 -(x(1)+1.0)^2)
  202.         exp1 = exp(-x(0)^2 -x(1)^2)
  203.         exp2 = exp(-(x(0)+1.0)^2 -x(1)^2)
  204.         xy   = x(0)/5.0 - x(0)^3 - x(1)^5
  205.         RETURN, DOUBLE( [-6.0*exp0*( x(0) * (1.0-x(0))^2 + (1.0-x(0)) ) $
  206.            -10.0 * exp1 * (-2.0 * x(0) * xy + (0.2 - 3.0 * x(0)^2))  $
  207.            +2.0 * (x(0) + 1.0) * exp2/3.0,                           $
  208.            -6.0 * (x(1) + 1.0) * (1.0 - x(0))^2 * exp0               $
  209.            -10.0 * exp1 * (-2.0 * x(1) * xy - (5.0 * x(1)^4))        $
  210.             +2.0 * x(1) * exp2/3.0] )
  211.     endif
  212.  
  213.     ;  Function # 3 .
  214.     ;
  215.     if (surfaceindex eq 3) then begin
  216.         exp0 = DOUBLE( exp(-x(0)^2 - x(1)^2) )
  217.         RETURN, DOUBLE([exp0 * (1.0 - 2.0*x(0)*x(0) - 4.0*x(0)*x(1) ),  $  
  218.             exp0 * (2.0 - 4.0*x(1)*x(1) - 2.0*x(0)*x(1) ) ] )
  219.     endif
  220.  
  221. end                   ; Of GRAD
  222.  
  223. ;---------------------------------------------------------------------------
  224. ;
  225. ;  Purpose: Execure the line search.
  226. ;
  227. pro LNSRCH, $
  228.     n, $                         ; IN: Dimension of the function
  229.     xold, $                      ; IN: Old point coordinates
  230.     fold, $                      ; IN: Old function value at x old
  231.     g, $                         ; IN/OUT: Vector of lenght n+1, 
  232.                                  ;         where g(0) is the step limit,  
  233.                                  ;         9(1) to g(n) is the returned
  234.                                  ;         gradient at point x
  235.     p,  $                        ; OUT: Newton direction (dim = n+1)
  236.     x, $                         ; OUT: New x  coordinates (dim = n+1),
  237.                                  ;      x(0) contains the check flag.
  238.     f, $                         ; OUT: Function value at point x
  239.     Surfaceindex = surfaceindex  ; IN ; function index
  240.  
  241.     ;  Make variable as double types.
  242.     ;
  243.     xold = double(xold)
  244.     fold = double(fold)
  245.     g = double(g)
  246.  
  247.     ;  Initialize constants.
  248.     ; 
  249.     stpmax = g(0)   ; Maximum number of steps
  250.     check = x(0)    ; check flag is false (0) on normal exit
  251.     tolx = 1.0D-7   ; convergence criterion on DELTA x
  252.     alf = 1.0D-4    ; Ensures sufficient decrease in function x
  253.  
  254.     ;  Check for proper dimension on x.
  255.     ;
  256.     sizen = SIZE(xold)
  257.     ns = sizen(1)
  258.     if (ns NE (n+1)) then begin
  259.         Message,'Error in LNSRCH, dimension of xold  NE n'
  260.     ENDIF
  261.  
  262.     check = 0.0D0
  263.     x(0) = check
  264.     p(0) = 0.0D0
  265.     sum = 0.0D0
  266.     sum = TOTAL(DOUBLE(p(1:n)*p(1:n)) )
  267.     sum = SQRT(sum)
  268.  
  269.     ;  Scale in attempted step is too big.
  270.     ;
  271.     if (sum GT stpmax) then begin
  272.         p(1:n) = DOUBLE(p(1:n)) * DOUBLE( stpmax)/sum
  273.     endif
  274.  
  275.     slope = 0.0D0
  276.     slope = TOTAL(double( g(1:n)*p(1:n) ))
  277.     test = 0.0D0
  278.  
  279.     for i = 1,n do begin
  280.         temp = ABS(p(i)) / MAX( [ABS(xold(i)), 1.0D0] )
  281.         if (temp GT test) THEN test = temp
  282.     endfor
  283.  
  284.     alamin = DOUBLE( tolx/test )
  285.     alam = 1.0D0
  286.     alam2 =0.0D0
  287.     fold2 =0.0D0
  288.     f2 = 0.0D0
  289.     count =0
  290.     dummy = 50
  291.  
  292.     ;  Start of iteration loop.
  293.     ;
  294.     while (dummy EQ 50 ) do begin
  295.         count = count + 1
  296.         x(1:n) = DOUBLE( xold(1:n) + alam*p(1:n) )
  297.         f = DOUBLE(FUNCOP(x(1:n),surfaceindex))
  298.  
  299.         if (alam LT alamin) then begin  ; if 1
  300.             x(1:n) = DOUBLE( xold(1:n) )
  301.             check = 1.0D0
  302.             x(0) = 1.0D0
  303.             RETURN
  304.         endif else if (f LE (fold+alf*alam*slope)) then begin ; if 1
  305.             RETURN
  306.         endif else begin     ; if 1 still
  307.             if ((alam EQ 1.0D0) ) then begin ; if 2
  308.                 tmplam = (-1.0D0*slope)/(2.0D0*(f-fold-slope))
  309.             endif else begin ; if 2
  310.                 rhs1 = DOUBLE( f -fold-alam*slope )
  311.                 rhs2 = DOUBLE( f2-fold2-alam2*slope )
  312.                 alamp2 = DOUBLE( alam*alam )
  313.                 alam2p2 = DOUBLE( alam2 * alam2 )
  314.                 alamdi = DOUBLE( alam-alam2 )
  315.                 a = DOUBLE( ((rhs1/alamp2)-(rhs2/alam2p2))/alamdi )
  316.                 b = DOUBLE( ( (-1.0D0*alam2*rhs1/alamp2) + $
  317.                     (alam*rhs2/alam2p2))/alamdi )
  318.  
  319.                 if (a EQ 0.0D0) then begin ; if 3
  320.                     tmplam =  -1.0D0*slope/(2.0D0*b)
  321.                 endif else begin ; if 3
  322.                     disc = b*b-3.0D0*a*slope
  323.  
  324.                     if (disc LT 0.0D0) then begin ; if 4
  325.                         Message, ' Roundoff problem in lnsrch'
  326.                     endif else begin
  327.                         tmplam = (-1.0D0*b + sqrt(disc))/(3.0D0*a)
  328.                     endelse   ; end of if loop 4
  329.                 endelse   ; end of if loop 3
  330.  
  331.                 if (tmplam GT 0.5D0*alam) THEN tmplam = 0.5D0*alam
  332.             endelse   ; end of if loop 2
  333.         endelse   ; end of if loop 1
  334.  
  335.         alam2 = DOUBLE(alam)
  336.         f2 = DOUBLE( f )
  337.         fold2 = DOUBLE( fold )
  338.         alam = MAX( [tmplam, 0.1D0*alam])
  339.     endwhile
  340.  
  341. end       ; of lnsrch procedure
  342.  
  343. ;---------------------------------------------------------------------------
  344. ;
  345. ;  Purpose:  This function computes the minimum point of a function
  346. ;            using the Dorovan-fletcher-Powell algorithm.
  347. ;            The IDL equivalent is called DFPMIN. The difference is
  348. ;            that DFP returns the intermediary steps.
  349. ;            For details and comments see Numerical recipes in C 2nd ed.
  350. ;
  351. pro dfp, $
  352.     n, $                         ; IN: Dimension of the function
  353.     p, $                         : IN: Vector (n+1) containing the initial 
  354.     gtol, $                      : IN: Limit of tolerance
  355.     iter, $                      : OUT: Number of iterations
  356.     xtot, $                      : OUT: Array ( iter+1, n+1) 
  357.                                  ;      that contains the data path.
  358.     Surfaceindex                 ; IN ; function index
  359.    
  360.     ;  Initialize constants.
  361.     ;
  362.     itmax = 200               ; maximum number of iterations
  363.     eps = 3.0D-7              ; machine precision
  364.     tolx = double(4.0* eps)   ; Convergence criterion on x
  365.     stpmx = double(100.0)     ; scaled max. step length allowed in LNSRCH
  366.  
  367.     ;  Initialize working arrays.
  368.     ;
  369.     dg = DBLARR(n+1)
  370.     g = DBLARR(n+1)
  371.     pnew = DBLARR(n+1)
  372.     hdg = DBLARR(n+1)
  373.     hessin = DBLARR(n+1,n+1)
  374.     xi = DBLARR(n+1)
  375.  
  376.     ;  Evaluate the function value at the starting point.
  377.     ;
  378.     fp = FUNCOP(double(p(1:n)), surfaceindex)
  379.  
  380.     ;  Evaluate the function gradient at the starting point.
  381.     ;  
  382.     g(1:n) = GRAD(DOUBLE(p(1:n)), surfaceindex)
  383.  
  384.     xtot(0,0:1) =DOUBLE( p(1:2) )
  385.     xtot(0,2) = DOUBLE(fp)
  386.  
  387.     ; Initialize the matrix hessin to the identity matrix
  388.     ;
  389.     hessin = hessin + DOUBLE(0.0)
  390.     for i = 1, n do begin
  391.         hessin(i,i) = DOUBLE(1.0)
  392.     endfor
  393.  
  394.     sum = TOTAL(DOUBLE( p(1:n)*p(1:n) ) )
  395.     xi(1:n) =  DOUBLE( -g(1:n) )
  396.     stpmax = stpmx*MAX([DOUBLE(n),SQRT(sum)])
  397.  
  398.     ; Begin loop over max number of iterations
  399.     ;
  400.     for its=1, itmax do begin
  401.         g(0) = DOUBLE(stpmax)
  402.         xi(0) = DOUBLE(0.0)
  403.         iter = its
  404.  
  405.         ;  Performs a line search by calling the function LNSRCH.
  406.         ;
  407.         LNSRCH,n,p,fp,g,xi,pnew,fret, Surfaceindex = surfaceindex
  408.    
  409.         xtot(its,0:1) = DOUBLE( pnew(1:2) )
  410.         xtot(its,2) = DOUBLE(fret )
  411.         check = DOUBLE( pnew(0) )
  412.         fp = DOUBLE(fret)
  413.         xi(1:n) = DOUBLE( pnew(1:n) - p(1:n) )
  414.         p(1:n) = pnew(1:n)
  415.         test = DOUBLE(0.0)
  416.         
  417.         for i = 1, n do begin
  418.            temp = DOUBLE(ABS(xi(i)) / MAX( [ABS(p(i)),double(1.0)]))
  419.            if (temp GT test) then test = temp
  420.         endfor
  421.    
  422.         if (test LT tolx) then RETURN
  423.  
  424.         dg(1:n) = DOUBLE( g(1:n) )
  425.    
  426.         ;  Get the dradient at the new point.
  427.         ;
  428.         g(1:n) = GRAD(DOUBLE( p(1:n)),surfaceindex)
  429.    
  430.         test = DOUBLE(0.0)
  431.         den = MAX([DOUBLE(fret),1.0D0])
  432.  
  433.         for i = 1, n do begin
  434.            temp = ABS(g(i))*MAX([ABS(p(i)),1.0D0])/den
  435.            IF(temp GT test) THEN test = temp
  436.         endfor
  437.  
  438.         ;  Return if test is smaller that the tolerance limit.
  439.         ;
  440.         if (test LT gtol) then RETURN
  441.  
  442.         dg(1:n) = double(  g(1:n) - dg(1:n) )
  443.         hdg(1:n) = double( dg(1:n) - dg(1:n) )
  444.         hdg(1:n) = double( dg(1:n) ## hessin(1:n,1:n) )
  445.  
  446.         ;  Initialize working variables.
  447.         ;
  448.         fac   = 0.0D0
  449.         fae   = 0.0D0
  450.         sumdg = 0.0D0
  451.         sumxi = 0.0D0
  452.  
  453.         fac   = TOTAL( DOUBLE( dg(1:n) *xi(1:n) ))
  454.         fae   = TOTAL(DOUBLE( dg(1:n) * hdg(1:n) ))
  455.         sumdg = TOTAL( DOUBLE( dg(1:n)*dg(1:n) ))
  456.         sumxi = TOTAL(DOUBLE( xi(1:n)*xi(1:n) ))
  457.  
  458.         if ((fac * fac) GT (eps*sumdg*sumxi) ) then begin
  459.            fac = 1.0D0/fac
  460.            fad = 1.0D0/fae
  461.            dg(1:n) = fac*xi(1:n) - fad*hdg(1:n)
  462.            hessin(1:n,1:n) = hessin(1:n,1:n) + $
  463.            transpose(fac*xi(1:n)##xi(1:n) - $
  464.            fad*hdg(1:n)##hdg(1:n) +fae*dg(1:n)##dg(1:n) )
  465.         endif
  466.  
  467.         xi = double( xi-xi )
  468.         hessint = double( -1.0D0*hessin )
  469.         xi(1:n) = double( g(1:n) ## hessint(1:n,1:n) )
  470.  
  471.     ;  End of iteration loop.
  472.     ;
  473.     endfor
  474.  
  475.     ; If the function does not returns from inside the
  476.     ; previous loop, then returns the message
  477.     ; 
  478.     MESSAGE, ' Too many iterations in DFPMIN'
  479.  
  480. end                  ; Of DFP function
  481.  
  482. ;
  483. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  484. ;       Unauthorized reproduction prohibited.
  485. ;
  486. ;+
  487. ;  FILE:
  488. ;       poweop.pro
  489. ;
  490. ;  CALLING SEQUENCE: poweop
  491. ;                    Called by : linmop.pro
  492. ;
  493. ;  PURPOSE:  Computes the minimum (optimization) of a function
  494. ;            using the POWELL algorthm.
  495. ;            Specific use for optim5.pro
  496. ;
  497. ;  MAJOR TOPICS: Data analysis
  498. ;
  499. ;  CATEGORY:
  500. ;       IDL 4.0
  501. ;
  502. ;  INTERNAL FUNCTIONS and PROCEDURES:
  503. ;       NONE.
  504. ;
  505. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  506. ;       pro funcop              - Returns the function value
  507. ;       pro linmop              - returns minimum along a direction.
  508. ;
  509. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  510. ;             Numerical recipes in C, 2nd ed.
  511. ;             Press, Vetterling, Teukolsky, and Flannery
  512. ;             Cambridge University Press
  513. ;             ISBN 0-521-43108-5
  514. ;
  515. ;  NAMED STRUCTURES:
  516. ;       NONE.
  517. ;
  518. ;  COMMON BLOCS:
  519. ;       NONE.
  520. ;
  521. ;  MODIFICATION HISTORY: Written by DAT,RSI, June 1996
  522. ;-
  523. ;---------------------------------------------------------------------------
  524. ;
  525. ;  Purpose: Compute the minimum of a function
  526. ;           using hte POWELL algortihm.
  527. ;
  528. pro POWEOP, $
  529.     p, $                     ; IN/OUT: initial starting point
  530.     xi, $                    ; IN: initial matrix (usually identity matrix)
  531.     ftol, $                  ; IN: fractional tolerance
  532.     iter, $                  ; OUT: number of iterations
  533.     fret, $                  ; OUT; function value at minimum point (fret)
  534.     path, $                  ; OUT: array containing the iteration steps.
  535.     status, $                : OUT: function status (0=failure, 1=success)
  536.     SURFACEINDEX=surfaceindex; IN: Surface index (1, 2,or 3)
  537.  
  538.     ;  Get the dimension of the function.
  539.     ;
  540.     n = p(0)
  541.     n = 2
  542.  
  543.     ;  Initialize working variables and arrays.
  544.     ;
  545.     ITMAX    = 200                  ; maximuim number of iterations allowed
  546.     NMAX     = 21                   ; maximum dimension of the system.
  547.     i        = 0                    ; integer
  548.     ibig     = 0                    ; integer
  549.     j        = 0                    ; integer
  550.     del      = 0.0D0                ; DOUBLE precision
  551.     fp       = 0.0D0                ; DOUBLE precision
  552.     fptt     = 0.0D0                ; DOUBLE precision
  553.     t        = 0.0D0                ; DOUBLE precision
  554.     pt       = DBLARR(nmax)         ; array
  555.     ptt      = DBLARR(nmax)         ; array
  556.     xit      = DBLARR(nmax)         ; array
  557.     status   = 1                    ; Status is Valid (OK)
  558.     iter     = 1                    ; iteration number
  559.     path     = DBLARR(itmax+1, n+1) ; array
  560.  
  561.  
  562.     ;  Get the function value at point p.
  563.     ;
  564.     fret = DOUBLE(FUNCOP(p(1:n), surfaceIndex) )
  565.  
  566.     ;  Initialize the vector path.
  567.     ;
  568.     path(0, 0:n-1 ) = p(1:n)
  569.     path(0, n) = fret
  570.  
  571.     ;  Set pt to p. Save the initial point.
  572.     ;
  573.     pt(1:n) = DOUBLE( p(1:n) )
  574.  
  575.     ;  Loop over iteration.
  576.     ;  Integer iter has already been set to one. It is
  577.     ;  incremented at the end of the for loop.
  578.     ;
  579.     dummy = 0
  580.  
  581.     for iter = 1, ITMAX+1 do begin
  582.  
  583.         fp = DOUBLE( fret )
  584.         ibig = 0
  585.         del = 0.0D0       ; Will be the biggest function decrease
  586.  
  587.         ;  In each iteration, loop over all directions in the set.
  588.         ;
  589.         for i = 1, n  do begin
  590.  
  591.             ;  Copy the direction.
  592.             ; 
  593.             xit(1:n) = DOUBLE( xi( 1:n, i)  )
  594.             fptt = DOUBLE(  fret )
  595.  
  596.             ;  Minimize along that direction.
  597.             ;
  598.             linmop, p, xit, n, fret, check, $
  599.                 SURFACEINDEX=surfaceIndex
  600.   
  601.             ;  Verify the check status of linmin.
  602.             ;
  603.             if (check EQ 0) then begin
  604.                 PRINT, ' Error in linmop NO 1 in POWEOP.pro'
  605.                 RETURN
  606.             endif
  607.  
  608.             ;  Augment the vector path.
  609.             ;
  610.             path(iter, 0:n-1) =DOUBLE( p(1:n) )
  611.             path(iter, n) = DOUBLE( fret )
  612.  
  613.             ;  Record it and if it is the largest decrease so far.
  614.             ;
  615.             if (ABS(fptt - fret) GT del ) then begin
  616.                 del = DOUBLE( ABS( fptt - fret) )
  617.                 ibig = i
  618.             endif
  619.  
  620.         endfor
  621.  
  622.         ;  Termination criterion.
  623.         ;
  624.         if ((2.0* (ABS(fp - fret)) ) LE $
  625.             (ftol*(ABS(fp) + ABS(fret)) ) ) then begin
  626.             status = 1  ; Valid
  627.             RETURN
  628.         endif
  629.  
  630.         if (iter  EQ ITMAX) then begin
  631.             status = 0   ; NOT Valid
  632.             PRINT, 'Number of iteration exceeds the limit in POWELM'
  633.             RETURN
  634.         endif
  635.  
  636.         ;  Construct the extrapolated point and the average
  637.         ;  direction moved. Save the old starting point.
  638.         ;
  639.         ptt(1:n) = DOUBLE( 2.0 * p(1:n) - pt(1:n) )
  640.         xit(1:n) = DOUBLE( p(1:n) - pt(1:n) )
  641.         pt(1:n) = DOUBLE( p(1:n) )
  642.  
  643.         ;  Function value at extrapolated point.
  644.         ;
  645.         fptt = DOUBLE( FUNCOP( ptt(1:n), surfaceIndex) )
  646.  
  647.  
  648.         ;  Move to the minimum of the new direction,
  649.         ;  and save the new direction.
  650.         ;
  651.         if (fptt LT fp) then begin
  652.             square1 = DOUBLE(  (fp - fret - del) * (fp - fret - del) )
  653.             square2 = DOUBLE( (fp - fptt) * (fp - fptt) )
  654.             t = DOUBLE( 2.0*(fp - (2.0*fret) + fptt) * square1 - $ 
  655.                 del * square2 )
  656.             if (t LT 0.0d) then begin
  657.                 linmop, p, xit, n, fret, check, $
  658.                     SURFACEINDEX=surfaceIndex
  659.  
  660.                 ;  Verify the check status of linmin.
  661.                 ;
  662.                 if (check EQ 0) then begin
  663.                     PRINT, ' Error in linmop NO 2 in POWEOP.pro'
  664.                     RETURN
  665.                 endif
  666.  
  667.                 ;  Augment or reset  the vector path.
  668.                 ;
  669.                 path(iter, 0:n-1) = p(1:n)
  670.                 path(iter, n) = DOUBLE( fret )
  671.                 xi(1:n, ibig) = DOUBLE( xi( 1:n, n) )
  672.                 xi(1:n, n)    = DOUBLE( xit(1:n) )
  673.             endif
  674.         endif
  675.     endfor
  676.  
  677. end   ; of POWELM.pro
  678.  
  679.  
  680. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  681. ;       Unauthorized reproduction prohibited.
  682. ;
  683. ;+
  684. ;  FILE:
  685. ;       sign.pro
  686. ;
  687. ;  CALLING SEQUENCE: sign
  688. ;                    Called by : linmop.pro
  689. ;                                breop5.pro
  690. ;
  691. ;  PURPOSE:  This function isolates the minimum to a fractional
  692. ;            precision of about tol using Brent's method.
  693. ;            Specific use for optim5.pro
  694. ;
  695. ;  MAJOR TOPICS: Data analysis
  696. ;
  697. ;  CATEGORY:
  698. ;       IDL 4.0
  699. ;
  700. ;  INTERNAL FUNCTIONS and PROCEDURES:
  701. ;        pro sign               - main procedure
  702. ;
  703. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  704. ;      NONE
  705. ;
  706. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  707. ;             Numerical recipes in C, 2nd ed.
  708. ;             Press, Vetterling, Teukolsky, and Flannery
  709. ;             Cambridge University Press
  710. ;             ISBN 0-521-43108-5
  711. ;
  712. ;  NAMED STRUCTURES:
  713. ;       none.
  714. ;
  715. ;  COMMON BLOCS:
  716. ;       none.
  717. ;
  718. ;  MODIFICATION HISTORY:Written by DAT,RSI, June 1996
  719. ;-
  720. ;---------------------------------------------------------------------------
  721. ;
  722. ;  Purpose:  This function returns  
  723. ;            absolute of arg1         if arg2 is positive
  724. ;            minus  absolute of arg1  if arg2 is less or equal to 0
  725. ;
  726. function SIGN, $
  727.     arg1, $       ; IN: argument 1
  728.     arg2, $       ; IN: argument2
  729.     STATUS=status ; OUT: 0 = failure, 1= success
  730.  
  731.     if (N_Elements(status) ) then begin
  732.         status = 1
  733.     endif
  734.  
  735.     nparam = N_Params()
  736.     if (nparam NE 2) THEN BEGIN
  737.         print, ' Error calling SIGN, must have 2 arguments'
  738.         status = 0
  739.         RETURN, double( arg1)
  740.     endif
  741.  
  742.     if (arg2 GT 0.0) then begin 
  743.         RETURN, double( ABS(arg1) )
  744.     endif else begin
  745.         RETURN, double( -(ABS(arg1))  )
  746.     endelse
  747.  
  748. end   ; of function SIGN
  749.  
  750.  
  751.  
  752. ;
  753. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  754. ;       Unauthorized reproduction prohibited.
  755. ;
  756. ;+
  757. ;  FILE:
  758. ;       f1dop5.pro
  759. ;
  760. ;  CALLING SEQUENCE: f1dop5
  761. ;                    Called by : breop5.pro
  762. ;
  763. ;  PURPOSE:  Returns the function value at the updated point.
  764. ;            Specific use for optim5.pro
  765. ;
  766. ;  MAJOR TOPICS: Data analysis
  767. ;
  768. ;  CATEGORY:
  769. ;       IDL 4.0
  770. ;
  771. ;  INTERNAL FUNCTIONS and PROCEDURES:
  772. ;       pro f1dop5              - Main procedure
  773. ;
  774. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  775. ;       pro funcop              - Returns the function value
  776. ;
  777. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  778. ;             Numerical recipes in C, 2nd ed.
  779. ;             Press, Vetterling, Teukolsky, and Flannery
  780. ;             Cambridge University Press
  781. ;             ISBN 0-521-43108-5
  782. ;
  783. ;  NAMED STRUCTURES:
  784. ;       none.
  785. ;
  786. ;  COMMON BLOCS:
  787. ;       none.
  788. ;
  789. ;  MODIFICATION HISTORY:Written by DAT,RSI, June 1996
  790. ;-
  791. ;---------------------------------------------------------------------------
  792. ;
  793. ;  Purpose: Returns the function value at the updated point.
  794. ;
  795. function f1dop5, $
  796.     x , $                       ; IN: non updated point 
  797.     pcom, $                     ; IN: vector for updaing the point
  798.     xicom, $                    ; IN: vector for updaing the point
  799.     ncom, $                     ; IN: dimension of the function
  800.     status, $                   ; OUT: status check flag.
  801.     SURFACEINDEX=surfaceindex   ; IN: surface index
  802.  
  803.     ;  Check the validity of input parameters.
  804.     ;
  805.     nparam = N_Params()
  806.     if (nparam NE 5) then begin
  807.         PRINT, ' Error calling f1dim, must have 5 arguments'
  808.         status = 0
  809.         RETURN, x
  810.     endif
  811.  
  812.     ;  Set status to valid.
  813.     ;
  814.     status = 1
  815.  
  816.     ;  Initialize variables.
  817.     ;
  818.     f  =  0.0D0            ; double precision
  819.     xt =  dblarr(ncom+1)   ; double precision
  820.  
  821.     ;  Update the point.
  822.     ;
  823.     xt(1:ncom) = double( pcom(1:ncom) + x*xicom(1:ncom) )
  824.  
  825.     ;  Get the function value at point  xt.
  826.     ;
  827.     f = double( FUNCOP( xt(1:ncom),surfaceindex ) )
  828.  
  829.     RETURN, f
  830.  
  831. end   ; of f1dim.pro
  832.  
  833.  
  834.  
  835. ;
  836. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  837. ;       Unauthorized reproduction prohibited.
  838. ;
  839. ;+
  840. ;  FILE:
  841. ;       breop5.pro
  842. ;
  843. ;  CALLING SEQUENCE: breop5
  844. ;                    Called by : linmop.pro
  845. ;
  846. ;  PURPOSE:  This function isolates the minimum to a fractional
  847. ;            precision of about tol using Brent's method.
  848. ;            Specific use for optim5.pro
  849. ;
  850. ;  MAJOR TOPICS: Data analysis
  851. ;
  852. ;  CATEGORY:
  853. ;       IDL 4.0
  854. ;
  855. ;  INTERNAL FUNCTIONS and PROCEDURES:
  856. ;       fun f1dop5              - Compute line vlaue.
  857. ;       pro sign                - Return the sign.
  858. ;       pro breop5              - Main procedure
  859. ;
  860. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  861. ;       pro funcop              - Returns the function value
  862. ;
  863. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  864. ;             Numerical recipes in C, 2nd ed.
  865. ;             Press, Vetterling, Teukolsky, and Flannery
  866. ;             Cambridge University Press
  867. ;             ISBN 0-521-43108-5
  868. ;
  869. ;  NAMED STRUCTURES:
  870. ;       none.
  871. ;
  872. ;  COMMON BLOCS:
  873. ;       none.
  874. ;
  875. ;  MODIFICATION HISTORY:Written by DAT,RSI, June 1996
  876. ;-
  877. ;---------------------------------------------------------------------------
  878. ;
  879. ;  Purpose: Isolate the minimum value.
  880. ;
  881. function BREOP5, $
  882.     ax, $                       ; IN: Inital abscissa
  883.     bx, $                       ; IN: Inital abscissa
  884.     cx, $                       ; IN: Inital abscissa
  885.     tol, $                      ; IN: Limit of tolerance
  886.     xmin, $                     ; OUT: Returned minimum 
  887.     status ,$                   ; OUT: 0 = Not valid for function BRENT
  888.     Ppassed = pp, $             ; IN: Vector passed to BRENT
  889.     Xipassed = xip, $           ; IN: Vector passed to BRENT
  890.     Npassed = np, $             ; IN: Value  passed to BRENT
  891.     SurfaceIndex=surfaceindex
  892.  
  893.     ;  Verify that p, n, and xi have been passed correctly.
  894.     ;
  895.     keywordp = N_Elements(pp)
  896.     if (keywordp EQ 0 ) then begin
  897.         PRINT, ' Error : p not passed to BRENT'
  898.         status = 0
  899.         RETURN, 0
  900.     endif
  901.  
  902.     keywordn = N_Elements(np)
  903.     if (keywordn EQ 0 ) then begin
  904.         PRINT, ' Error : n not passed to BRENT'
  905.         status = 0
  906.         RETURN, 0
  907.     endif
  908.  
  909.     keywordxi = N_Elements(xip)
  910.     if (keywordxi EQ 0 ) then begin
  911.         PRINT, ' Error : xi not passed to BRENT'
  912.         status = 0
  913.         RETURN, 0
  914.     endif
  915.  
  916.     ;  Initialize working variables and arrays.
  917.     ;
  918.     itmax  = 100           ; maximum number of iterations
  919.     cgold  = 0.381966D0   ; golden ratio 
  920.     zeps   = 1.0D-10       ; small number for numerical protection
  921.     iter   = 0             ; iteration number
  922.     a      = 0.0D0         ; DOUBLE precision
  923.     b      = 0.0D0         ; DOUBLE precision
  924.     d      = 0.0D0         ; DOUBLE precision
  925.     etemp  = 0.0D0         ; DOUBLE precision
  926.     fu     = 0.0D0         ; DOUBLE precision
  927.     fv     = 0.0D0         ; DOUBLE precision
  928.     fw     = 0.0D0         ; DOUBLE precision
  929.     fx     = 0.0D0         ; DOUBLE precision
  930.     p      = 0.0D0         ; DOUBLE precision
  931.     q      = 0.0D0         ; DOUBLE precision
  932.     r      = 0.0D0         ; DOUBLE precision
  933.     tol1   = 0.0D0         ; DOUBLE precision
  934.     tol2   = 0.0D0         ; DOUBLE precision
  935.     u      = 0.0D0         ; DOUBLE precision
  936.     v      = 0.0D0         ; DOUBLE precision
  937.     w      = 0.0D0         ; DOUBLE precision
  938.     x      = 0.0D0         ; DOUBLE precision
  939.     xm     = 0.0D0         ; DOUBLE precision
  940.     e      = 0.0D0         ; Distance moved on the step before last
  941.     status = 1             ; Status is Valid (OK)
  942.  
  943.     if (ax LT cx) then a = ax else a = cx
  944.     if (ax GT cx) then b = ax else b = cx
  945.  
  946.     x = DOUBLE(bx)
  947.     w = DOUBLE(bx)
  948.     v = DOUBLE(bx)
  949.  
  950.     ;  Call the function f1dim.pro  .
  951.     ;
  952.     fx =DOUBLE( f1dop5(x, pp, xip, np , $
  953.         checkF1dim, SURFACEINDEX=surfaceindex) )
  954.  
  955.     ;  Verify the check status of f1dop5.
  956.     ;
  957.     if (checkF1dim EQ 0) then begin
  958.         PRINT, ' error in F1dim function NO 1 (in BRENT)'
  959.         status = 0
  960.         RETURN, 0
  961.     endif
  962.  
  963.     fv = DOUBLE(fx)
  964.     fw = DOUBLE(fx)
  965.  
  966.     ;  Main loop program.
  967.     ;
  968.     for iter = 1, ITMAX do begin
  969.         xm =  DOUBLE( 0.5*(a + b) )
  970.         tol1 = DOUBLE( tol*ABS(x) + zeps )
  971.         tol2 = DOUBLE( 2.0 * tol1 ) 
  972.    
  973.         ;  Test for done here.
  974.         ;
  975.         if (ABS( x - xm) LE ( tol2 - 0.5d*(b-a)) ) then begin
  976.            xmin = x
  977.            RETURN, fx
  978.         endif
  979.  
  980.         ;  Construct a trial parabolic fit.
  981.         ;
  982.         if (ABS(e) GT tol1 ) then begin
  983.             r = DOUBLE( (x - w) * (fx - fv) )
  984.             q =  DOUBLE( (x - v) * (fx - fw) )
  985.             p =  DOUBLE( (x - v)*q - (x - w)*r )
  986.             q =  DOUBLE( 2.0d * (q - r) )
  987.             if (q GT 0.0d) then p =  DOUBLE( -p )
  988.             q = DOUBLE( ABS(q) )
  989.             etemp = DOUBLE( e )
  990.             e =  DOUBLE( d )
  991.      
  992.             ;  Check for the acceptability for the parabolic fit
  993.             ;  otherwise, take the golden section step into the larger
  994.             ;  of the 2 segments.
  995.             if (( ABS(p) GE ABS(0.5d*q*etemp) ) OR $
  996.                 ( p LE q*(a - x) ) OR $
  997.                 ( p GE q*(b - x) ) ) then begin
  998.  
  999.                 if (x GE xm ) then begin
  1000.                     e = DOUBLE(a - x)
  1001.                 endif else begin
  1002.                     e = DOUBLE(b - x)
  1003.                 endelse
  1004.           
  1005.                 d = DOUBLE( cgold * e )
  1006.    
  1007.             endif else begin
  1008.  
  1009.                 d = DOUBLE( p/q )
  1010.                 u = DOUBLE( x + d )
  1011.                 if (( (u-a) LT tol2 ) OR ( (b - u) LT tol2 ) ) then begin
  1012.                     d = DOUBLE( SIGN(tol1, (xm - x), Status = status ) )
  1013.                     if (status EQ 0) then begin
  1014.                         print, ' error in SIGN function NO 1 (in BRENT)'
  1015.                         RETURN, 0
  1016.                     endif
  1017.                 endif
  1018.             endelse
  1019.         endif else begin
  1020.  
  1021.             if (x GE xm ) then begin
  1022.                 e = DOUBLE(a - x)
  1023.             endif else begin
  1024.                 e = DOUBLE(b - x)
  1025.             endelse
  1026.  
  1027.             d = DOUBLE( cgold * e )
  1028.         endelse
  1029.  
  1030.         if (ABS(d) GE tol1) then begin
  1031.             u =  DOUBLE( x+d)
  1032.         endif else begin
  1033.             u = x + SIGN( tol1, d, Status = status)
  1034.             if (status EQ 0) then begin
  1035.                print, ' error in SIGN function NO 2 (in BRENT)'
  1036.                RETURN, 0
  1037.             endif
  1038.         endelse
  1039.    
  1040.         ;  Call the function f1dop5.pro.
  1041.         ;
  1042.         fu = f1dop5(u, pp, xip, np , $
  1043.             checkF1dim, SURFACEINDEX=surfaceindex)
  1044.  
  1045.         ;  Verify the check status of f1dop5.
  1046.         ;
  1047.         if (checkF1dim EQ 0) then begin
  1048.             PRINT, ' error in F1dop5 function NO 2 (in BRENT)'
  1049.             status = 0
  1050.             RETURN, 0
  1051.         endif
  1052.  
  1053.         ;  This is the one function evaluation per iteration.
  1054.         ;
  1055.         if ( fu LE fx) then begin
  1056.             if (u GE x ) then  a = DOUBLE(x)  else b=DOUBLE(x)
  1057.             v = w
  1058.             w = x
  1059.             x = u
  1060.             fv = fw
  1061.             fw = fx
  1062.             fx = fu
  1063.         endif else begin
  1064.  
  1065.             if (u LT x) then a=u else b=u
  1066.  
  1067.             if ((fu LE fw) OR (w EQ x) ) then begin
  1068.                 v  =  DOUBLE(w)
  1069.                 w  = DOUBLE(u)
  1070.                 fv = DOUBLE(fw)
  1071.                 fw = DOUBLE(fu)
  1072.             endif else if ( (fu LE fv) OR (v EQ x) OR (v EQ w) )then begin
  1073.                 v  = DOUBLE(u)
  1074.                 fv = DOUBLE(fu)
  1075.             endif
  1076.         endelse
  1077.  
  1078.     endfor  ; end of the main loop
  1079.  
  1080.     ;  If the program reaches here,
  1081.     ;  then the number of iteration is too high.
  1082.     ;
  1083.     status = 0
  1084.     PRINT, ' Error in BRENT, too many iterations
  1085.     RETURN, 0
  1086.  
  1087. end ; of BRENT
  1088.  
  1089. ;
  1090. ;  Copyright (c) 1997, Research Systems, Inc. All rights reserved.
  1091. ;       Unauthorized reproduction prohibited.
  1092. ;
  1093. ;+
  1094. ;  FILE:
  1095. ;       mnbop5.pro
  1096. ;
  1097. ;  CALLING SEQUENCE: mnbop5
  1098. ;                    Called by : linmop.pro
  1099. ;
  1100. ;  PURPOSE:  This function isolates the minimum to a fractional
  1101. ;            precision of about tol using Brent's method.
  1102. ;            Specific use for optim5.pro
  1103. ;
  1104. ;  MAJOR TOPICS: Data analysis
  1105. ;
  1106. ;  CATEGORY:
  1107. ;       IDL 4.0
  1108. ;
  1109. ;  INTERNAL FUNCTIONS and PROCEDURES:
  1110. ;       fun f1dop5              - Compute line value.
  1111. ;       pro sign                - Return the sign.
  1112. ;       pro breop5              - Main procedure
  1113. ;
  1114. ;  EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
  1115. ;       pro funcop              - Returns the function value
  1116. ;
  1117. ;  REFERENCE: IDL Reference Guide, IDL User's Guide
  1118. ;             Numerical recipes in C, 2nd ed.
  1119. ;             Press, Vetterling, Teukolsky, and Flannery
  1120. ;             Cambridge University Press
  1121. ;             ISBN 0-521-43108-5
  1122. ;
  1123. ;  NAMED STRUCTURES:
  1124. ;       none.
  1125. ;
  1126. ;  COMMON BLOCS:
  1127. ;       none.
  1128. ;
  1129. ;  MODIFICATION HISTORY:Written by DAT,RSI, June 1996
  1130. ;-
  1131. ;---------------------------------------------------------------------------
  1132. ;
  1133. ;  Purpose:  This function searches in the downhill direction,  
  1134. ;            and returns the new points ax, bx, cx (as parametrs) that bracket 
  1135. ;            a minimum of the function f1dop5. The status is returned
  1136. ;            as value by the function.
  1137. ;
  1138. function MNBOP5, $
  1139.     ax, $                       ; IN/OUT: initial/final point 
  1140.     bx, $                       ; IN/OUT: initial/final point
  1141.     cx, $                       ; OUT: final point
  1142.     fa, $                       ; OUT: value function at ax
  1143.     fb, $                       ; OUT: value function at bx
  1144.     fc, $                       ; OUT: value function at cx
  1145.     PPASSED = pp, $             ; IN: parameter used by f1dop5
  1146.     XIPASSED = xip, $           ; IN: parameter used by f1dop5
  1147.     NPASSED = np, $             ; IN: parameter used by f1dop5
  1148.     SURFACEINDEX = surfaceIndex ; IN: surface index (1, 2, or 3)
  1149.  
  1150.     ;  Verify that pp, np, and xip have been passed. 
  1151.     ;
  1152.     keywordp = N_ELEMENTS(pp)
  1153.     if (keywordp EQ 0 ) then begin
  1154.         PRINT, ' Error in MNBOP5 : p not passed to MNBOP5'
  1155.         status = 0
  1156.         RETURN, status
  1157.     endif
  1158.  
  1159.     keywordn = N_ELEMENTS(np)
  1160.     if (keywordn EQ 0 ) then begin
  1161.         PRINT, ' Error in MNBOP5 : n not passed to MNBOP5'
  1162.         status = 0
  1163.         RETURN, 0
  1164.     endif
  1165.  
  1166.     keywordxi = N_ELEMENTS(xip)
  1167.     if (keywordxi EQ 0 ) then begin
  1168.         PRINT, ' Error in MNBOP5 : xi not passed to MNBOP5'
  1169.         status = 0
  1170.         RETURN, 0
  1171.     endif
  1172.  
  1173.     ;  Initialize working variables and arrays.
  1174.     ;
  1175.     gold   = 1.618034d     ; default ratio 
  1176.     glimit = 100.0D0       ; max. magnification for parabolic-fit step
  1177.     tiny   = 1.0D-20       ; double precision
  1178.     ulim   = 0.0D0         ; double precision
  1179.     u      = 0.0D0         ; double precision
  1180.     r      = 0.0D0         ; double precision
  1181.     q      = 0.0D0         ; double precision
  1182.     fu     = 0.0D0         ; double precision
  1183.     dum    = 0.0D0         ; double precision
  1184.     status = 1             ; Status is Valid (OK)
  1185.  
  1186.     ;  Call the function f1dop5.pro twice .
  1187.     ;
  1188.     fa = f1dop5(ax, pp, xip, np , checkF1dim, SURFACEINDEX=surfaceindex)
  1189.  
  1190.     ;  Verify the check status of f1dop5.
  1191.     ;
  1192.     if (checkF1dim EQ 0) then begin
  1193.         PRINT, ' error in F1dop5 function NO 1 (in MNBOP5)'
  1194.         status = 0
  1195.         RETURN, 0
  1196.     endif
  1197.  
  1198.     fb = f1dop5(bx, pp, xip, np , checkF1dim, SURFACEINDEX=surfaceindex)
  1199.  
  1200.     ;  Verify the check status of f1dim.
  1201.     ;
  1202.     if (checkF1dim EQ 0) then begin
  1203.         PRINT, ' error in F1dop5 function NO 2 (in MNBOP5)'
  1204.         status = 0
  1205.         RETURN, 0
  1206.     endif
  1207.  
  1208.  
  1209.     ;  Switch roles of a and b so that we can go
  1210.     ;  downhill in the direction for a to b.
  1211.     ;
  1212.     if (fb GT fa) then begin
  1213.         dum = ax
  1214.         ax = bx
  1215.         bx = dum
  1216.         dum = fb
  1217.         fb = fa
  1218.         fa = dum
  1219.     endif
  1220.  
  1221.  
  1222.     ;  First guess for c.
  1223.     ;
  1224.     cx = DOUBLE( bx + gold*(bx - ax) )
  1225.  
  1226.     fc = f1dop5(cx, pp, xip, np , checkF1dim, SURFACEINDEX=surfaceindex)
  1227.  
  1228.     ;  Verify the check status of f1dop5.
  1229.     ;
  1230.     if (checkF1dim EQ 0) then begin
  1231.         PRINT, ' error in F1dim function NO 3 (in MNBOP5)'
  1232.         status = 0
  1233.         RETURN, 0
  1234.     endif
  1235.  
  1236.  
  1237.     ;  Keep returning here until we bracket.
  1238.     ;
  1239.     while (fb GT fc) do begin
  1240.  
  1241.         ;  Compute u by parabolic extrapolation from a,b,c.
  1242.         ;  Tiny is used to prevent any possible division by zero.
  1243.         ;
  1244.         r = DOUBLE( (bx - ax) * (fb - fc) )
  1245.         q = DOUBLE( (bx - cx) * (fb - fa) )
  1246.         valueFmax =  DOUBLE( MAX( ABS(q-r), tiny) )
  1247.         valueSign = SIGN( valueFmax,  (q-r), STATUS=status)
  1248.         if (status eq 0) then begin
  1249.             PRINT, ' Error calling SIGN (in MNBOP5)'
  1250.             RETURN, 0
  1251.         endif
  1252.  
  1253.         u = DOUBLE( bx - ( (bx-cx)*q - (bx-ax)*r ) / $
  1254.             (2.0 * valueSign ) ) 
  1255.  
  1256.         ulim =  DOUBLE( bx + glimit*(cx - bx) )
  1257.    
  1258.         ;  Test various possibility.
  1259.         ;
  1260.         if (((bx - u)*(u - cx)) GT 0.0 )  then begin
  1261.  
  1262.             fu = f1dop5(u, pp, xip, np ,checkF1dim, $
  1263.                 SURFACEINDEX=surfaceindex)
  1264.  
  1265.             ;  Verify the check status of f1dim.
  1266.             ;
  1267.             if (checkF1dim EQ 0) then begin
  1268.                PRINT, ' error in F1dop5 function NO 4 (in MNBOP5)'
  1269.                status = 0
  1270.                RETURN, 0
  1271.             endif
  1272.  
  1273.             ;  Got a minimum between b and c.
  1274.             ;
  1275.             if (fu LT fc) then begin
  1276.                 ax =  DOUBLE( bx )
  1277.                 bx =  DOUBLE( u )
  1278.                 fa =  DOUBLE( fb )
  1279.                 fb =  DOUBLE( fu )
  1280.                 RETURN, 1   ; status should be 1 (ok) here
  1281.  
  1282.             ;  Got a minimum between a and u.
  1283.             ;
  1284.             endif else if (fu GT fb) then begin
  1285.                 cx = DOUBLE( u )
  1286.                 fc = DOUBLE( fu )
  1287.                 RETURN, 1   ; status should be 1 (ok) here
  1288.             endif
  1289.  
  1290.             ;  Parabolic fit was no use, Use default magnification.
  1291.             ;
  1292.             u =  double( cx + gold*(cx - bx) )
  1293.  
  1294.             fu = f1dop5(u, pp, xip, np ,checkF1dim, $
  1295.                 SURFACEINDEX=surfaceindex)
  1296.  
  1297.             ;  Verify the check status of f1dim.
  1298.             ;
  1299.             if (checkF1dim EQ 0) then begin
  1300.                 PRINT, ' error in F1dop5 function NO 5 (in MNBOP5)'
  1301.                 status = 0
  1302.                 RETURN, 0
  1303.             endif
  1304.       
  1305.         ; Parabolic fit is between c and its allowed limit.
  1306.         ;
  1307.         endif else if (((cx - u)*( u - ulim)) GT 0.0) then begin
  1308.  
  1309.             fu = f1dop5(u, pp, xip, np ,checkF1dim, $
  1310.                 SURFACEINDEX=surfaceindex)
  1311.  
  1312.             ;  Verify the check status of f1dim.
  1313.             ;
  1314.             if (checkF1dim EQ 0) then begin
  1315.                 PRINT, ' error in f1dop5 function NO 6 (in MNBOP5)'
  1316.                 status = 0
  1317.                 RETURN, 0
  1318.             endif
  1319.       
  1320.             if (fu LT fc) then begin
  1321.                 arg4 = cx+gold*(cx-bx)
  1322.                 bx = cx
  1323.                 cx = u
  1324.                 u = arg4
  1325.                 fut = f1dop5(u, pp, xip, np , checkF1dim, $
  1326.                     SURFACEINDEX=surfaceindex)
  1327.  
  1328.                 if (checkF1dim EQ 0) then begin
  1329.                     PRINT, ' error in f1dop5 function NO 7 (in MNBOP5)'
  1330.                     status = 0
  1331.                     RETURN, 0
  1332.                 endif
  1333.          
  1334.                 fb = fc
  1335.                 fc = fu
  1336.                 fu = fut
  1337.             endif
  1338.  
  1339.         ;  Limit parabolic u to maximum allowed limit.
  1340.         ;
  1341.         endif else if (( u - ulim)*(ulim - cx) GE 0.0) then begin
  1342.             u = DOUBLE( ulim )
  1343.  
  1344.             fu = f1dop5(u, pp, xip, np , checkF1dim, $
  1345.                 SURFACEINDEX=surfaceindex)
  1346.  
  1347.             ;  Verify the check status of f1dim.
  1348.             ;
  1349.             if (checkF1dim EQ 0) then begin
  1350.                 PRINT, ' error in F1dim function NO 8 (in MNBOP5)'
  1351.                 status = 0
  1352.                 RETURN, 0
  1353.              endif
  1354.  
  1355.         ;  Reject parabolic u, use default magnification.
  1356.         ;
  1357.         endif else begin
  1358.             u =  double( cx + gold*(cx - bx) )
  1359.  
  1360.             fu = f1dop5(u, pp, xip, np , checkF1dim, $
  1361.                 SURFACEINDEX=surfaceindex)
  1362.  
  1363.             ;  Verify the check status of f1dim.
  1364.             ;
  1365.             if (checkF1dim EQ 0) then begin
  1366.                 PRINT, ' error in F1dim function NO 9 (in MNBOP5)'
  1367.                 status = 0
  1368.                 RETURN, 0
  1369.              endif
  1370.         endelse
  1371.  
  1372.         ;  Eliminate oldest point and continue.
  1373.         ;
  1374.         ax = bx
  1375.         bx = cx
  1376.         cx = u
  1377.         fa = fb
  1378.         fb = fc
  1379.         fc = fu
  1380.     endwhile
  1381.  
  1382.     RETURN, status
  1383.  
  1384. end  ; of mnbrak.pro
  1385.       
  1386. pro LINMOP, p, xi, n , fret, status,SurfaceIndex = surfaceindex
  1387. ;
  1388. ; TAKEN from Numerical Recipes in C (2nd ed.).
  1389. ;
  1390. ; This procedure finds the minimum point along the 
  1391. ; direction xi and passing by the point p
  1392. ;
  1393. ; INPUT
  1394. ;      xi             :   inital direction ( dim : n, xi(1) to xi(n)  )
  1395. ;      n              :   The dimension of the surface minus
  1396. ;                          one. Example : the surface z = 3 x + 6y +2.5 
  1397. ;                         (which is a plane) has n = 2 .
  1398. ; surfaceindex :   index that indicates the selected surface (1,2,3)
  1399. ;
  1400. ; INPUT/OUTPUT
  1401. ;      p          :   Initial starting point stored in
  1402. ;                     p(1) to p(n).
  1403. ;                     The returned value is the
  1404. ;                     location of the function minimum along the direction
  1405. ;                     given by xi
  1406. ; OUTPUT
  1407. ;      fret         :   value of the function at point p (minimum)
  1408. ;      status     :   0 = Not valid 
  1409. ;
  1410. ;
  1411. ; THIS PROCEDURE IS CALLED BY
  1412. ;      Procedure POWELM.PRO
  1413. ;
  1414. ; THIS PROCEDURE CALLS  FOR
  1415. ;      Function BRENT.pro
  1416. ;      Function MNBRAK.pro
  1417. ;      Function F1dim.pro ( read this as f one dim. pro)
  1418. ;      Function FUNC.pro 
  1419.  
  1420.  
  1421.     ;  Initialize working variables and arrays.
  1422.     ;
  1423.     tol    = 2.0D-6        ; tolerance passed to BRENT.pro
  1424.     xx     = 1.0D0         ; double precision
  1425.     xmin   = 0.0D0         ; double precision
  1426.     fx     = 0.0D0         ; double precision
  1427.     fb     = 0.0D0         ; double precision
  1428.     fa     = 0.0D0         ; double precision
  1429.     bx     = 0.0D0         ; double precision
  1430.     ax     = 0.0D0         ; double precision
  1431.     status = 1             ; Status is Valid (OK)
  1432.  
  1433.     psaved1 = double(p)
  1434.     psaved2 = double( p )
  1435.     xisaved1 = double( xi )
  1436.     xisaved2 = double( xi )
  1437.     nsaved1 = double( n)
  1438.     nsaved2 = double( n )
  1439.  
  1440.     ;  Call the function MNBRAK.pro .
  1441.     ;  IMPORTANT NOTE : p,xi and n MUST be passed,
  1442.     ;  otherwise an error will be issued.
  1443.     ;
  1444.     checkMnbrak = MNBOP5(ax, xx, bx, fa, fx, fb, $
  1445.         Ppassed = psaved1, Xipassed = xisaved1, Npassed = nsaved1, $
  1446.         SurfaceIndex = surfaceindex) 
  1447.  
  1448.     ;  Verify the check status of MNBRAK.
  1449.     ;
  1450.     if (checkMnbrak EQ 0) then begin
  1451.         print, ' error in MNBRAK NO 1 in LINMIN.pro'
  1452.         status = 0
  1453.         RETURN
  1454.     endif
  1455.  
  1456.     ;  Call the function BRENT.pro.
  1457.     ;  IMPORTANT NOTE : p,xi and n MUST be passed,
  1458.     ;  Otherwise an error will be issued.
  1459.     ;
  1460.     fret = double( BREOP5(ax, xx, bx, tol, xmin, checkBrent, $ 
  1461.         Ppassed = psaved2, Xipassed = xisaved2, Npassed = nsaved2, $
  1462.         SurfaceIndex = surfaceindex) )
  1463.  
  1464.     ;  Verify the check status of  BRENT.PRO.
  1465.     ;
  1466.     if (checkBrent EQ 0) then begin
  1467.         print, ' error in BRENT NO 1 in LINMIN.pro'
  1468.         status = 0
  1469.         RETURN
  1470.     endif
  1471.  
  1472.     ;  Construct the vector results to return.
  1473.     ;
  1474.     xi(1:n) = double( xi(1:n)* xmin )
  1475.     p(1:n) = double( p(1:n) + xi(1:n)) 
  1476.  
  1477. end  ; of LINMIN
  1478.  
  1479.  
  1480.  
  1481. ;---------------------------------------------------------------------------
  1482. ;
  1483. ;  Purpose:   Initialize 16 working (predefined) colors.
  1484. ;
  1485. pro opt_Init16colors, $
  1486.     colornames  ; IN: color anmes array
  1487.  
  1488.     TVLCT, 255, 255, 255, colornames(0)     ; white
  1489.     TVLCT, 255, 255,   0, colornames(1)     ; yellow
  1490.     TVLCT, 200, 180, 255, colornames(2)     ; lavender
  1491.     TVLCT,  80, 255, 255, colornames(3)     ; aqua
  1492.     TVLCT, 255, 100, 150, colornames(4)     ; pink
  1493.     TVLCT,  55, 255,  55, colornames(5)     ; green
  1494.     TVLCT, 255,  30,  30, colornames(6)     ; red
  1495.     TVLCT, 255, 180,  80, colornames(7)     ; orange
  1496.     TVLCT,  80,  80, 255, colornames(8)     ; blue
  1497.     TVLCT, 180, 180, 180, colornames(9)     ; lt_gray
  1498.     TVLCT,  55, 130, 100, colornames(10)    ; med_green
  1499.     TVLCT,  80,  30,  30, colornames(11)    ; brown
  1500.     TVLCT,  55,  55,   0, colornames(12)    ; olive
  1501.     TVLCT, 100,  30,  80, colornames(13)    ; purple
  1502.     TVLCT,  55,  55,  55, colornames(14)    ; dk_gray
  1503.     TVLCT,   0,   0,   0, colornames(15)    ; black
  1504.  
  1505. end  ; of opt_init16colors
  1506.  
  1507. ; -----------------------------------------------------------------------------
  1508. ;
  1509. ;  Purpose:  This procedure initialize the starting points
  1510. ;            (4 points per surface, 3 surfaces, total 12 points).
  1511. ;
  1512. pro opt_InitPoint, $
  1513.     startpoints   ; OUT: start points locations.
  1514.  
  1515.     point = dblarr(2)        ; working array
  1516.     i = 1                    ; surface index
  1517.  
  1518.     point  = [ -1.8d, 2.0d ]
  1519.     startpoints(i,1, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1520.  
  1521.     point  = [ -1.0d, 1.9d ]
  1522.     startpoints(i,2, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1523.  
  1524.     point  = [ 0.0d, 1.9d ]
  1525.     startpoints(i,3, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1526.  
  1527.     point  = [ 1.0d, 2.0d ]
  1528.     startpoints(i,4, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1529.  
  1530.     i = 2
  1531.  
  1532.     point  = [ -1.6d, 1.0d ]
  1533.     startpoints(i,1, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1534.  
  1535.     point  = [ -0.7d, 0.7d ]
  1536.     startpoints(i,2, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1537.  
  1538.     point  = [ 0.75d, -0.875d ]
  1539.     startpoints(i,3, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1540.  
  1541.     point  = [ 1.375d, -1.8d ]
  1542.     startpoints(i,4, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1543.  
  1544.     i = 3
  1545.  
  1546.     point  = [ -0.9d, 0.7d ]
  1547.     startpoints(i,1, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1548.  
  1549.     point  = [ -0.2d, 0.5d ]
  1550.     startpoints(i,2, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1551.  
  1552.     point  = [ 0.2d, -0.2d ]
  1553.     startpoints(i,3, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1554.  
  1555.     point  = [ 0.6d, -0.5d ]
  1556.     startpoints(i,4, 1:3) = [point(0), point(1), FUNCOP(point(0:1),i) ]
  1557.  
  1558. end  ; of INITPOINTS
  1559.  
  1560. ; -----------------------------------------------------------------------------
  1561. ;
  1562. ;  Purpose:  This procedure resets the starting points.
  1563. ;
  1564. pro opt_ResetPoint, $
  1565.     event,  $     ; IN: event structure.
  1566.     pointIndex    ; IN: point index
  1567.  
  1568.     ;  Get the info structure.
  1569.     ;
  1570.     WIDGET_CONTROL, event.top, GET_UVALUE=info, /NO_COPY
  1571.  
  1572.     info.pointindex = pointindex
  1573.  
  1574.     ;  Reset the results labels.
  1575.     ;
  1576.     opt_ResetResult, info.ResultLBLx,info.ResultLBLy, $
  1577.         info.ResultLBLz,info.ResultLBLi
  1578.  
  1579.     if (info.surfaceIndex EQ 1) then begin
  1580.         saved3Dt = info.saved3d1
  1581.         pix = info.pixmap1ID
  1582.     endif else if (info.surfaceIndex EQ 2) then begin
  1583.         saved3Dt = info.saved3d2
  1584.         pix = info.pixmap2ID
  1585.     endif else begin
  1586.         saved3Dt = info.saved3d3
  1587.         pix = info.pixmap3ID
  1588.     endelse
  1589.  
  1590.     ;  Desensitize all the points buttons.
  1591.     ;
  1592.     WIDGET_CONTROL, info.Point1ID, SENSITIVE=0
  1593.     WIDGET_CONTROL, info.Point2ID, SENSITIVE=0
  1594.     WIDGET_CONTROL, info.Point3ID, SENSITIVE=0
  1595.     WIDGET_CONTROL, info.Point4ID, SENSITIVE=0
  1596.  
  1597.     ;  Redraw the surface from the pixmap.
  1598.     ;
  1599.     WSET, info.drawWinID
  1600.     Device, Copy=[ 0,0,info.drawxsize, $
  1601.         info.drawysize, 0,0, pix]
  1602.  
  1603.     ;  Plot all the points.
  1604.     ;
  1605.     opt_PlotPoint, info.surfaceIndex, pointindex, $
  1606.         info.charscale,info.red, $
  1607.         saved3dt, info.startPoints,  $
  1608.         info.startPointsdev, info.fontflag
  1609.  
  1610.     ;  Sensitize the appropriate buttons
  1611.     ; and change the information text label.
  1612.     ;
  1613.     WIDGET_CONTROL, info.startbuttonID, SENSITIVE=1
  1614.     WIDGET_CONTROL, info.surface1ID, SENSITIVE=1
  1615.     WIDGET_CONTROL, info.surface2ID, SENSITIVE=1
  1616.     WIDGET_CONTROL, info.surface3ID, SENSITIVE=1
  1617.     WIDGET_CONTROL, info.StartButtonID, SENSITIVE=1
  1618.     WIDGET_CONTROL, info.FindButton, SENSITIVE=1
  1619.  
  1620.     if (pointindex EQ 1) then begin
  1621.         WIDGET_CONTROL, info.Point1ID, SENSITIVE=0
  1622.         WIDGET_CONTROL, info.Point2ID, SENSITIVE=1
  1623.         WIDGET_CONTROL, info.Point3ID, SENSITIVE=1
  1624.         WIDGET_CONTROL, info.Point4ID, SENSITIVE=1
  1625.     endif else if (pointindex EQ 2) then begin
  1626.         WIDGET_CONTROL, info.Point1ID, SENSITIVE=1
  1627.         WIDGET_CONTROL, info.Point2ID, SENSITIVE=0
  1628.         WIDGET_CONTROL, info.Point3ID, SENSITIVE=1
  1629.         WIDGET_CONTROL, info.Point4ID, SENSITIVE=1
  1630.     endif else if (pointindex EQ 3) then begin
  1631.         WIDGET_CONTROL, info.Point1ID, SENSITIVE=1
  1632.         WIDGET_CONTROL, info.Point2ID, SENSITIVE=1
  1633.         WIDGET_CONTROL, info.Point3ID, SENSITIVE=0
  1634.         WIDGET_CONTROL, info.Point4ID, SENSITIVE=1
  1635.     endif else begin
  1636.         WIDGET_CONTROL, info.Point1ID, SENSITIVE=1
  1637.         WIDGET_CONTROL, info.Point2ID, SENSITIVE=1
  1638.         WIDGET_CONTROL, info.Point3ID, SENSITIVE=1
  1639.         WIDGET_CONTROL, info.Point4ID, SENSITIVE=0
  1640.     endelse
  1641.  
  1642.     ;  Add the annotation string on top of the points.
  1643.     ;
  1644.     opt_ShowStart, info.surfaceIndex, pointindex, $
  1645.         info.StartPoints, info.StartLBL
  1646.  
  1647.     ;  Restore the info structure.
  1648.     ;
  1649.     WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1650.  
  1651. end  ; opt_ResetPoint
  1652.  
  1653. ; -----------------------------------------------------------------------------
  1654. ;
  1655. ;  Purpose:   Redraw the appropriate surface.
  1656. ;
  1657. pro opt_ResetSurface, $
  1658.     event, $       ; IN: event structure
  1659.     surfaceindex   ; IN: surface index (1,2, or 3)
  1660.  
  1661.     ;  Get the info structure.
  1662.     ;
  1663.     WIDGET_CONTROL, event.top, GET_UVALUE=info, /NO_COPY
  1664.  
  1665.     ;  Desensitize all the surface buttons.
  1666.     ;
  1667.     WIDGET_CONTROL, info.surface1ID, sensitive = 0
  1668.     WIDGET_CONTROL, info.surface2ID, sensitive = 0
  1669.     WIDGET_CONTROL, info.surface3ID, sensitive = 0
  1670.  
  1671.     ;  Reset the results labels.
  1672.     ;
  1673.     opt_ResetResult, info.ResultLBLx,info.ResultLBLy, $
  1674.         info.ResultLBLz,info.ResultLBLi
  1675.  
  1676.     info.surfaceIndex = surfaceindex
  1677.     info.pointIndex = 0
  1678.     WSET, info.drawWinID
  1679.  
  1680.     ;  Reset the data system coordinates to the appropriate
  1681.     ;  surface by creating an empty plot....
  1682.     ;
  1683.     if (surfaceindex EQ 1) then begin
  1684.         pix = info.pixmap1ID
  1685.         saved3dt = info.saved3d1
  1686.         SURFACE, info.zdata1, info.xdata1, $
  1687.             info.ydata1,/nodata,  $
  1688.             AZ=10, AX=60, $
  1689.             XSTYLE=4, YSTYLE=4, ZSTYLE=4, $
  1690.             POSITION=[0.1, 0.20, 0.95, 1.]
  1691.     endif else if (surfaceindex EQ 2) then begin
  1692.         pix = info.pixmap2ID
  1693.         saved3dt = info.saved3d2
  1694.         SURFACE, info.zdata2, info.xdata2, $
  1695.             info.ydata2,  /nodata, $
  1696.             AZ=30, AX=30,$
  1697.             XRANGE=[-2,2], $
  1698.             YRANGE=[-2,2], $
  1699.             ZRANGE=[-7,7],$
  1700.             XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  1701.             POSITION=[0.05, 0.05, 0.95, 1.]
  1702.     endif else begin
  1703.         pix = info.pixmap3ID
  1704.         saved3dt = info.saved3d3
  1705.         SURFACE, info.zdata3, info.xdata3,$
  1706.             info.ydata3, /SAVE,/NODATA, $
  1707.             AZ=30,AX= 30, $
  1708.             XRANGE=[-1,1], $
  1709.             YRANGE=[-1,1], $
  1710.             ZRANGE=[-1,1],$
  1711.             XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  1712.             POSITION=[0.05, 0.05, 0.95, 1.]
  1713.     endelse
  1714.  
  1715.     ;  Copy the selected surface from 
  1716.     ;  the pixmap into the drawing window.
  1717.     ;
  1718.     DEVICE, COPY=[0, 0, info.drawxsize, $
  1719.         info.drawysize, 0, 0, pix]
  1720.  
  1721.     ;  Plot the starting points.
  1722.     ;
  1723.     opt_PlotPoint, surfaceindex, 1, info.charscale,info.red, $
  1724.         saved3dt, info.startPoints, $
  1725.         info.startPointsdev, info.fontflag,  /ALL
  1726.  
  1727.     ;  Sensitize the appropriate buttons
  1728.     ;  and reset the starting point and information labels.
  1729.     ;
  1730.     WIDGET_CONTROL, info.startbuttonID, SENSITIVE=1
  1731.     WIDGET_CONTROL, info.StartButtonID, SENSITIVE=1
  1732.     WIDGET_CONTROL, info.Point1ID, SENSITIVE=1
  1733.     WIDGET_CONTROL, info.Point2ID, SENSITIVE=1
  1734.     WIDGET_CONTROL, info.Point3ID, SENSITIVE=1
  1735.     WIDGET_CONTROL, info.Point4ID, SENSITIVE=1
  1736.     WIDGET_CONTROL, info.FindButton, SENSITIVE=0
  1737.     WIDGET_CONTROL, info.StartLBL, $
  1738.         Set_Value = 'X = ----   Y = ----'
  1739.  
  1740.     strings = '  Select another surface or a starting point'
  1741.  
  1742.     if (surfaceindex EQ 1) then begin
  1743.         WIDGET_CONTROL, info.surface2ID, SENSITIVE=1
  1744.         WIDGET_CONTROL, info.surface3ID, SENSITIVE=1
  1745.     endif else if( surfaceindex EQ 2) then begin
  1746.         WIDGET_CONTROL, info.surface1ID, SENSITIVE=1
  1747.         WIDGET_CONTROL, info.surface3ID, SENSITIVE=1
  1748.     endif else if( surfaceindex EQ 3) then begin
  1749.         WIDGET_CONTROL, info.surface1ID, SENSITIVE=1
  1750.         WIDGET_CONTROL, info.surface2ID, SENSITIVE=1
  1751.     endif
  1752.  
  1753.     ;  Restore the info structure.
  1754.     ;
  1755.     WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1756.  
  1757. end ; of RESETSURFACE
  1758.  
  1759.  
  1760. ; -----------------------------------------------------------------------------
  1761. ;
  1762. ;  Purpose:   This procedure open a text window and display
  1763. ;             the text file.
  1764. ;
  1765. pro About_Optimization, $
  1766.     group_leader    ;  IN: group leader identifer
  1767.  
  1768.     ;  Return if the info window has already been opened
  1769.     ;
  1770.     if (Xregistered('XDisplayFile') NE 0)  then return
  1771.  
  1772.     XDisplayFile, filepath("optimize.txt", $
  1773.         SUBDIR=['examples','demo','demotext']), $
  1774.         DONE_BUTTON='Done', $
  1775.         TITLE="About Optimization", $
  1776.         GROUP=group_leader, WIDTH=55, HEIGHT=14
  1777.  
  1778. end   ; of About_Optimization
  1779.  
  1780. ; -----------------------------------------------------------------------------
  1781. ;
  1782. ;  Purpose:   Reset the result labels.
  1783. ;
  1784. pro opt_ResetResult, $
  1785.     LabeLxID, $    ; IN: x label
  1786.     LabeLyID, $    ; IN: y label
  1787.     LabeLzID, $    ; IN: z label
  1788.     LabeLiID       ; IN: iteration label
  1789.  
  1790.     ; Reset the result labels
  1791.     ;
  1792.     WIDGET_CONTROL, LabeLxID, $
  1793.         SET_VALUE= 'X =        -----    -----'
  1794.     WIDGET_CONTROL, LabeLyID, $
  1795.         SET_VALUE= 'Y =        -----    -----'
  1796.     WIDGET_CONTROL, LabeLzID, $
  1797.         SET_VALUE= 'Z =        -----    -----'
  1798.     WIDGET_CONTROL, LabeLiID, $
  1799.         SET_VALUE= '#STEPS =   --       --'
  1800.  
  1801. end  ; of opt_ResetResult
  1802.  
  1803. ; -----------------------------------------------------------------------------
  1804. ;
  1805. ;  Purpose: Shows the value of the selected starting point in
  1806. ;           the starting point label.
  1807. ;
  1808. PRO opt_ShowStart, $
  1809.     surfaceIndex, $   ; IN: surface index 
  1810.     pointIndex, $     ; IN: point index
  1811.     startPoints, $    ; IN: starting points array
  1812.     StartLBL          ; IN: start label identifier
  1813.  
  1814.     ;  Shows the value of the selected starting point in
  1815.     ;  the starting point label.
  1816.     ;
  1817.     xval = StartPoints(surfaceIndex, pointIndex, 1)
  1818.     yval = StartPoints(surfaceIndex, pointIndex, 2)
  1819.     xs = STRING(xval, FORMAT='(f6.1)')
  1820.     xs = STRTRIM(xs,2)
  1821.     ys = STRING(yval, FORMAT='(f6.1)')
  1822.     ys = STRTRIM(ys,2)
  1823.     stringh ='X = ' + xs +'   Y = ' + ys
  1824.     WIDGET_CONTROL, StartLBL, SET_VALUE=stringh
  1825.  
  1826. END  ; opt_ShowStart
  1827.  
  1828. ; -----------------------------------------------------------------------------
  1829. ;
  1830. ;  Purpose:  This porcedure plots the starting points
  1831. ;            either the  selected one  or all of them.
  1832. ;
  1833. pro opt_PlotPoint, $
  1834.     surfaceIndex, $     ; IN: surface index
  1835.     pointnumber, $      ; IN: point index
  1836.     charscale, $        ; IN: character scaling factor
  1837.     red, $              ; IN: red color index in color table
  1838.     saved3dIndex, $     ; IN: 3-D transformation matrix
  1839.     startPoints, $      ; IN: start points (data coordinates)
  1840.     startPointsdev, $   ; IN: start points (device coordinates)
  1841.     fontFlag, $         ; IN: font( 0 = hardware, 1 = user selection)
  1842.     All=all             ; IN: (opt) plot all the point
  1843.  
  1844.     !P.T = 0
  1845.     si = surfaceIndex
  1846.  
  1847.     ; If all of the points
  1848.     ;
  1849.     if (N_Elements(all) EQ 1) then begin
  1850.         init =1
  1851.         final = 4
  1852.     endif else begin
  1853.         init = pointNumber
  1854.         final = pointNumber
  1855.     endelse
  1856.  
  1857.  
  1858.     for i = init, final do begin
  1859.         if (i EQ 1) then begin
  1860.             if (fontFlag EQ 0 ) then begin
  1861.                 string ='!51!x'
  1862.             endif else begin
  1863.                 string = '1'
  1864.             endelse
  1865.         endif else if (i EQ 2) then begin
  1866.             if (fontFlag EQ 0 ) then begin
  1867.                 string ='!52!x'
  1868.             endif else begin
  1869.                 string ='2'
  1870.             endelse
  1871.         endif else if (i EQ 3) then begin
  1872.             if (fontFlag EQ 0 ) then begin
  1873.                 string ='!53!x'
  1874.             endif else begin
  1875.                 string ='3'
  1876.             endelse
  1877.         endif else begin
  1878.             if (fontFlag EQ 0 ) then begin
  1879.                 string ='!54!x'
  1880.             endif else begin
  1881.                 string ='4'
  1882.             endelse
  1883.         endelse
  1884.  
  1885.         ;  Plot the point annotations (P1, P2, etc...).
  1886.         ;
  1887.         if (fontFlag EQ 1) then begin
  1888.             !P.FONT = 0
  1889.             DEVICE, FONT="TIMES*24"
  1890.         endif
  1891.  
  1892.         XYOUTS, startPointsdev(si,i,1), startPointsdev(si,i,2)+10, $
  1893.             string, COLOR=red, CHARSIZE=2.*charScale, $
  1894.             /DEVICE, ALIGNMENT=0.5
  1895.  
  1896.         if (fontFlag EQ 1) then begin
  1897.             DEVICE, FONT="TIMES*8"
  1898.             !P.FONT = -1
  1899.         endif
  1900.  
  1901.         ;  Set he 3-D transformation matrix.
  1902.         ;
  1903.         !P.T = saved3dIndex
  1904.  
  1905.         ;  Plot the point symbols.
  1906.         ;
  1907.         PLOTS, $
  1908.             [startPoints(si,i,1), startPoints(si,i,1), $
  1909.                 startPoints(si,i,1)], $
  1910.             [startPoints(si,i,2), startPoints(si,i,2), $
  1911.                 startPoints(si,i,2)], $
  1912.             [startPoints(si,i,3), startPoints(si,i,3), $
  1913.                 startPoints(si,i,3)], $
  1914.             PSYM=1, /Z, COLOR=red, THICK=2.6*charScale, /T3D
  1915.     endfor
  1916.  
  1917.  
  1918.     ;  Reset the 3d transformation.
  1919.     ;
  1920.     !P.T = 0
  1921.  
  1922. end  ; of opt_PlotPoint
  1923.  
  1924. ; -----------------------------------------------------------------------------
  1925. ;
  1926. ;  Purpose:  Main event handler.
  1927. ;
  1928. pro Optimize_event, $
  1929.     event      ; IN: event structure
  1930.  
  1931.     ;  Quit the application using the close box.
  1932.     ;
  1933.     if (TAG_NAMES(event, /STRUCTURE_NAME) EQ $
  1934.         'WIDGET_KILL_REQUEST') then begin
  1935.         WIDGET_CONTROL, event.top, /DESTROY
  1936.         RETURN
  1937.     endif
  1938.  
  1939.  
  1940.     ;  Get the info structure.
  1941.     ;
  1942.     WIDGET_CONTROL, event.top, GET_UVALUE=info, /NO_COPY
  1943.  
  1944.     ;  Get the button identifier value.
  1945.     ;
  1946.     WIDGET_CONTROL, event.id, GET_VALUE=buttonValue
  1947.  
  1948.     ;  Branch to the appropriate button.
  1949.     ;
  1950.     case buttonValue of
  1951.  
  1952.         'Surface 1' :  begin
  1953.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1954.             opt_ResetSurface, event, 1
  1955.         end  ; of  Surface 1
  1956.  
  1957.         'Surface 2' :  begin
  1958.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1959.             opt_ResetSurface, event, 2
  1960.         end  ; of  Surface 2
  1961.  
  1962.         'Surface 3' :  begin
  1963.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1964.             opt_ResetSurface, event, 3
  1965.         end  ; of  Surface 3
  1966.  
  1967.         'Reset Surface' : begin
  1968.             surfaceIndex = info.surfaceIndex
  1969.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1970.             opt_ResetSurface, event, surfaceIndex
  1971.         end  ; of  reset surface
  1972.  
  1973.         'Point 1' :  begin
  1974.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1975.             opt_ResetPoint, event, 1
  1976.         end  ; of Point 1
  1977.  
  1978.         'Point 2' :  begin
  1979.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1980.             opt_ResetPoint, event, 2
  1981.         end  ; of Point 2
  1982.  
  1983.         'Point 3' :  begin
  1984.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1985.             opt_ResetPoint, event, 3
  1986.         end  ; of Point 3
  1987.  
  1988.         'Point 4' :  begin
  1989.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1990.             opt_ResetPoint, event, 4
  1991.         end  ; of Point 4
  1992.  
  1993.         'Quit' :  begin
  1994.             WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  1995.             WIDGET_CONTROL, event.top, /DESTROY
  1996.         end   ;  of Quit
  1997.  
  1998.         'Find Minimum' : begin
  1999.    
  2000.              ;  Sensitize (desensitize) the appropeiate buttons
  2001.              ;  And reset the information label.
  2002.              ;
  2003.              WIDGET_CONTROL, info.StartButtonID, SENSITIVE=0
  2004.              WIDGET_CONTROL, info.SurfaceButtonID, SENSITIVE=0
  2005.              WIDGET_CONTROL, info.wResetButton, SENSITIVE=0
  2006.              WIDGET_CONTROL, info.findButton, SENSITIVE=0
  2007.              WIDGET_CONTROL, info.QuitButton, SENSITIVE=0
  2008.              WIDGET_CONTROL, info.AboutButtonID, SENSITIVE=0
  2009.  
  2010.              ;  Initialize variable for the DFP function call.
  2011.              ;
  2012.              itmax = 200                 ; Maximum number of iterations
  2013.              gtol = double(1.0D-7)       ; Limit of tolerance
  2014.              xtot = DBLARR(itmax+1,3)    ; Array that contains all the points
  2015.              n = 2                       ; Dimension of the banana function
  2016.              pin = dblarr(3)
  2017.              pin = [0.0, 0.0, 2.5]       ; Initial starting point (p(0) = 0.0)
  2018.              pin(1) = info.startpoints(info.surfaceIndex,$
  2019.                  info.pointIndex,1)
  2020.              pin(2) = info.startpoints(info.surfaceIndex, $
  2021.                  info.pointIndex,2)
  2022.              iter = 0
  2023.  
  2024.              ;  Computes the data points (path).
  2025.              ;
  2026.              DFP, n, pin, gtol, iter, xtot, info.surfaceIndex
  2027.  
  2028.              ;  Put the data path (contained in xtot) within 3 vectors.
  2029.              ;
  2030.              xvec = DBLARR(iter+1)
  2031.              yvec = DBLARR(iter+1)
  2032.              zvec = DBLARR(iter+1)
  2033.              xvec(0:iter) = xtot(0:iter,0)
  2034.              yvec(0:iter) = xtot(0:iter,1)
  2035.              zvec(0:iter) = xtot(0:iter,2)
  2036.  
  2037.              ;  Compute the path with the Powell algorithm.
  2038.              ;
  2039.              ;  Initialize variables for the Powell function call.
  2040.              ;
  2041.              itmax = 200                 ; Maximum number of iterations
  2042.              gtol = double(1.0D-7)       ; Limit of tolerance
  2043.              xtot = DBLARR(itmax+1,3)    ; Array that contains all the points
  2044.              n = 2                       ; Dimension of the banana function
  2045.              pin = dblarr(3)
  2046.              pin = [0.0, 0.0, 2.5]       ; Initial starting point (p(0) = 0.0)
  2047.              pin(0) = 2
  2048.              pin(1) = info.startpoints(info.surfaceIndex,$
  2049.                 info.pointIndex,1)
  2050.              pin(2) = info.startpoints(info.surfaceIndex, $
  2051.                 info.pointIndex,2)
  2052.              ftol = 1.0d-7
  2053.              xi = dblarr(3,3)
  2054.              xi = xi - xi
  2055.              xi(0,0) = 1.0d
  2056.              xi(1,1) = 1.0d
  2057.              xi(2,2) = 1.0d
  2058.  
  2059.              ;  Find the minimum using the POWELL algrithm.
  2060.              ;
  2061.              poweop, pin, xi, ftol, $
  2062.                 iterpo, fretpo,  pathpo, $
  2063.                 statuspo, Surfaceindex = info.surfaceindex
  2064.  
  2065.  
  2066.              ;  Put the data path (contained in pathpo) within 3 vectors.
  2067.              ;
  2068.              xvecpo = DBLARR(iterpo+1)
  2069.              yvecpo = DBLARR(iterpo+1)
  2070.              zvecpo = DBLARR(iterpo+1)
  2071.              xvecpo(0:iterpo) = pathpo(0:iterpo,0)
  2072.              yvecpo(0:iterpo) = pathpo(0:iterpo,1)
  2073.              zvecpo(0:iterpo) = pathpo(0:iterpo,2)
  2074.  
  2075.              ;  Now plotting the DFP only...
  2076.              ;
  2077.              WAIT, 1.0
  2078.              if (info.surfaceIndex EQ 1) then begin
  2079.                  !P.T = info.saved3d1
  2080.              endif else IF( info.surfaceIndex EQ 2) then begin
  2081.                  !P.T = info.saved3d2
  2082.              endif else  begin
  2083.                  !P.T = info.saved3d3
  2084.              endelse
  2085.  
  2086.              ;  Showing the message in the information label.
  2087.              ;
  2088.              string1 = '!3DFP!X'
  2089.              string2 = '!3Powell!X'
  2090.  
  2091.              if (info.fontflag EQ 0) then begin
  2092.                  string1 = '!3DFP!X'
  2093.                  string2 = '!3Powell!X'
  2094.              endif else if (info.fontflag EQ 1) then begin
  2095.                  !p.FONT = 0
  2096.                  Device, FONT = "TIMES*18"
  2097.                  string1 = 'DFP'
  2098.                  string2 = 'Powell'
  2099.              endif
  2100.  
  2101.              XYOUTS, 0.80, 0.1, $
  2102.                  string1, COLOR=info.orange, CHARSIZE=1.0*info.charScale, $
  2103.                  /NORMAL, CHARTHICK=1.0
  2104.  
  2105.              XYOUTS, 0.80, 0.05, $
  2106.                  string2, COLOR=info.white, CHARSIZE=1.0*info.charScale, $
  2107.                  /NORMAL, CHARTHICK=1.0
  2108.  
  2109.  
  2110.              if (info.fontflag EQ 1) then begin
  2111.                  DEVICE, FONT="TIMES*8"
  2112.                  !P.FONT = -1
  2113.              endif
  2114.  
  2115.              ;  Plot the path points, wait 0.3 seconds between each points.
  2116.              ;  Plot the DFP path first.
  2117.              ;
  2118.              for i = 0,iter-1 do begin
  2119.                  PLOTS, [xvec(i),xvec(i+1)], [yvec(i),yvec(i+1)], $
  2120.                      [zvec(i), zvec(i+1)], $
  2121.                      COLOR=info.orange, PSYM=-4, /T3D, $
  2122.                      THICK=2.25,SymSize = 1.6
  2123.                      WAIT, 0.3
  2124.              endfor
  2125.  
  2126.              WAIT, 1.0
  2127.  
  2128.              ;  Plot the POWELL path.
  2129.              ;
  2130.              for i = 0,iterpo-1 do begin
  2131.                  PLOTS, [xvecpo(i),xvecpo(i+1)], [yvecpo(i),yvecpo(i+1)], $
  2132.                      [zvecpo(i), zvecpo(i+1)], $
  2133.                      COLOR=info.white, PSYM=-4, /T3D, $
  2134.                      THICK=2.25,SymSize=1.6
  2135.                  WAIT, 0.3
  2136.              endfor
  2137.  
  2138.              ; Compute the device coordinates of the minimum point.
  2139.              ;
  2140.              mindev = DBLARR(3)
  2141.              mindev = $
  2142.                  CONVERT_COORD(xvec(iter-1), $
  2143.                  yvec(iter-1), $
  2144.                  zvec(iter-1), $
  2145.                  /to_device, /t3D)
  2146.  
  2147.              ;  Annotate the Minimum point.
  2148.              ;
  2149.              string = '!5Min!x'
  2150.              if (info.fontFlag EQ 1) then begin
  2151.                  !P.FONT = 0
  2152.                  DEVICE, FONT="TIMES*24"
  2153.                  string='Min'
  2154.              endif
  2155.              if (info.fontflag EQ 0) then begin
  2156.                  string = '!5Min!x'
  2157.              endif else if (info.fontflag EQ 1) then begin
  2158.                  !p.FONT = 0
  2159.                  DEVICE, FONT = "TIMES*24"
  2160.                  string = 'Min'
  2161.              endif
  2162.              XYOUTS,mindev(0)+ 20,mindev(1) -10,$
  2163.                  string, COLOR=info.red, $
  2164.                  CHARSIZE=2.*info.charScale,/device
  2165.              if (info.fontflag EQ 1) then begin
  2166.                  Device, FONT = "TIMES*8"
  2167.                  !p.FONT = -1
  2168.              endif
  2169.       
  2170.              ;  Draw a circle at the minimum point location.
  2171.              ;
  2172.              XYOUTS, xvec(iter-1), yvec(iter-1), Z=zvec(iter-1), $
  2173.                  '!20P!X', COLOR=info.red, CHARSIZE=2.*info.charScale, $
  2174.                  /T3D, CHARTHICK=2.0
  2175.  
  2176.              ;  Show the results in label box.
  2177.              ;
  2178.              xval = xvec(iter-1)
  2179.              yval = yvec(iter-1)
  2180.              zval = zvec(iter-1)
  2181.  
  2182.              xvalpo = xvecpo(iterpo-1)
  2183.              yvalpo = yvecpo(iterpo-1)
  2184.              zvalpo = zvecpo(iterpo-1)
  2185.  
  2186.              xsdfp = STRING(xval, FORMAT='(f7.4)')
  2187.              xsdfp = STRTRIM(xsdfp,2)
  2188.              ysdfp = STRING(yval, FORMAT='(f7.4)')
  2189.              ysdfp = STRTRIM(ysdfp,2)
  2190.              zsdfp = STRING(zval, FORMAT='(f7.4)')
  2191.              zsdfp = STRTRIM(zsdfp,2)
  2192.              isdfp = STRING(iter, FORMAT='(i3)')
  2193.              isdfp = STRTRIM(isdfp,2)
  2194.  
  2195.              xspow = STRING(xvalpo, FORMAT='(f7.4)')
  2196.              xspow = STRTRIM(xspow,2)
  2197.              yspow = STRING(yvalpo, FORMAT='(f7.4)')
  2198.              yspow = STRTRIM(yspow,2)
  2199.              zspow = STRING(zvalpo, FORMAT='(f7.4)')
  2200.              zspow = STRTRIM(zspow,2)
  2201.              ispow = STRING(iterpo, FORMAT='(i3)')
  2202.              ispow = STRTRIM(ispow,2)
  2203.  
  2204.              space = ' '
  2205.  
  2206.              if (xval   GE 0.0) then sp1 = space else sp1 =''
  2207.              if (xvalpo GE 0.0) then sp2 = space else sp2 =''
  2208.              stringx = 'X =     ' + sp1 + xsdfp + '    ' + sp2 + xspow
  2209.              if (yval   GE 0.0) then sp1 = space else sp1 =''
  2210.              if (yvalpo GE 0.0) then sp2 = space else sp2 =''
  2211.              stringy = 'Y =     ' + sp1 + ysdfp + '    ' + sp2 + yspow
  2212.              if (zval   GE 0.0) then sp1 = space else sp1 =''
  2213.              if (zvalpo GE 0.0) then sp2 = space else sp2 =''
  2214.              stringz = 'Z =     ' + sp1 + zsdfp + '    ' + sp2 + zspow
  2215.              stringi ='#STEPS = ' + isdfp + '            ' +ispow
  2216.  
  2217.  
  2218.              WIDGET_CONTROL, info.ResultLBLx, SET_VALUE=stringx
  2219.              WIDGET_CONTROL, info.ResultLBLy, SET_VALUE=stringy
  2220.              WIDGET_CONTROL, info.ResultLBLz, SET_VALUE=stringz
  2221.              WIDGET_CONTROL, info.ResultLBLi, SET_VALUE=stringi
  2222.  
  2223.              WIDGET_CONTROL, info.SurfaceButtonID, SENSITIVE=1
  2224.              WIDGET_CONTROL, info.StartButtonID, SENSITIVE=1
  2225.              WIDGET_CONTROL, info.findbutton, SENSITIVE=0
  2226.              WIDGET_CONTROL, info.wResetButton, SENSITIVE=1
  2227.              WIDGET_CONTROL, info.QuitButton, SENSITIVE=1
  2228.              WIDGET_CONTROL, info.AboutButtonID, SENSITIVE=1
  2229.  
  2230.              ;  Restore the info structure.
  2231.              ;
  2232.              WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  2233.  
  2234.          end  ; of Find Minimum
  2235.  
  2236.          'About Optimization' :  begin
  2237.  
  2238.              ;  Shows the text file into a window.
  2239.              ;
  2240.              About_Optimization, event.top
  2241.              WIDGET_CONTROL, event.top, SET_UVALUE=info, /NO_COPY
  2242.  
  2243.          end  ; of About optimization
  2244.  
  2245.     endcase   ;  of  buttonValue
  2246.  
  2247. end  ; Optimize_event
  2248.  
  2249. ; -----------------------------------------------------------------------------
  2250. ;
  2251. ;  Purpose:  Cleanup procedure.
  2252. ;
  2253. pro Optimize_Cleanup, $
  2254.     tlb     ;  IN: top level base
  2255.  
  2256.     ;  Get the info structure.
  2257.     ;
  2258.     WIDGET_CONTROL, tlb, GET_UVALUE=info, /no_copy
  2259.  
  2260.     ;  Restore the previous color table.
  2261.     ;
  2262.     TVLCT, info.colorTable
  2263.  
  2264.     ;  Delete the 3 pixmaps.
  2265.     ;
  2266.     WDELETE, info.pixmap1ID
  2267.     WDELETE, info.pixmap2ID
  2268.     WDELETE, info.pixmap3ID
  2269.  
  2270.     if WIDGET_INFO(info.groupBase, /VALID_ID) then $
  2271.         WIDGET_CONTROL, info.groupBase, /MAP
  2272.  
  2273. end  ; Of OPTIMIZE_Cleanup
  2274.  
  2275. ; -----------------------------------------------------------------------------
  2276. ;
  2277. ;  Purpose:  This application shows the minimization (optimization)
  2278. ;            methods by the DFP and POWELL algortihms.
  2279. ;            the user select a starting point and the path
  2280. ;            is displayed.
  2281. ;
  2282. pro D_Optimize, $
  2283.     GROUP=group, $     ; IN: (opt) group identifier
  2284.     APPTLB = appTLB    ; OUT: (opt) TLB of this application
  2285.  
  2286.     ;  This procedure plot the path taken by the
  2287.     ;  DFP and the Powell  minimzation algorithms.
  2288.     ;
  2289.     ;  The user can choose between 3 surfaces and 4 points
  2290.     ;  in each of these surfaces.
  2291.  
  2292.     ;  If there is an error, go back to the user.
  2293.     ;
  2294.     ON_Error, 1
  2295.  
  2296.  
  2297.     ;  Find the number of available colors.
  2298.     ;
  2299.     Window, /FREE, XSIZE = 2, YSIZE=2, /PIXMAP
  2300.     pixid = !D.WINDOW
  2301.     WDELETE, pixid
  2302.     numcolors = !D.TABLE_SIZE
  2303.  
  2304.     ;  If 'OPTIMIZE' is already open, then return...
  2305.     ;
  2306.     if (Xregistered('d_optimize') NE 0) then begin
  2307.       RETURN
  2308.     endif
  2309.  
  2310.     fontFlag = 0     ; 0 = use hersey fonts
  2311.                      ; 1 = use the buitlt in PS fonts (helvetica, etc..)
  2312.  
  2313.     ;  Save the current color table.
  2314.     ;
  2315.     TVLCT, savedR, savedG, savedB, /GET
  2316.     colorTable = [[savedR], [savedG], [savedB]]
  2317.  
  2318.     ; Check the validity of the group identifier.
  2319.     ;
  2320.     ngroup = N_ELEMENTS(group)
  2321.     if (ngroup NE 0) then begin
  2322.         check = WIDGET_INFO(group, /VALID_ID)
  2323.         if (check NE 1) then begin
  2324.             print,'Error, the group identifier is not valid'
  2325.             print, 'Return to the main application'
  2326.             RETURN
  2327.         endif
  2328.         groupBase = group
  2329.     endif else groupBase = 0L
  2330.  
  2331.     begintime = systime(1)
  2332.  
  2333.     ;  Get the screen size.
  2334.     ;
  2335.     DEVICE, GET_SCREEN_SIZE = screenSize
  2336.  
  2337.     ;  Create the starting up message.
  2338.     ;
  2339.     if (ngroup EQ 0) then begin
  2340.         drawbase = startmes()
  2341.     endif else begin
  2342.         drawbase = startmes(GROUP=group)
  2343.     endelse
  2344.  
  2345.     ;  Set up the color table.
  2346.     ;  Here we shift down the color table by 'nshift' positions.
  2347.     ;
  2348.     Window, /FREE, XSIZE=10, YSIZE=10, /PIXMAP
  2349.     pixid = !D.WINDOW
  2350.     WDELETE, pixid
  2351.     max_color = !D.TABLE_SIZE
  2352.  
  2353.     if (max_color GE 2L) then begin
  2354.         white = (max_color - 1L)  > 0L
  2355.         yellow = (max_color - 2L) > 0L
  2356.         lavender = (max_color - 3L) > 0L
  2357.         aqua = (max_color - 4L) > 0L
  2358.         pink = (max_color - 5L) > 0L
  2359.         green = (max_color - 6L) > 0L
  2360.         red = (max_color - 7L) > 0L
  2361.         orange = (max_color - 8L) > 0L
  2362.         blue = (max_color - 9L) > 0L
  2363.         lt_gray = (max_color - 10L) > 0L
  2364.         med_green = (max_color - 11L) > 0L
  2365.         brown = (max_color - 12L) > 0L
  2366.         olive = (max_color - 13L) > 0L
  2367.         purple = (max_color - 14L) > 0L
  2368.         dk_gray = (max_color - 15L) > 0L
  2369.         black = (max_color - 16L) > 0L
  2370.     endif else begin
  2371.         white = 1L
  2372.         yellow = 1L
  2373.         lavender = 1L
  2374.         aqua = 1L
  2375.         pink = 1L
  2376.         green = 1L
  2377.         red = 1L
  2378.         orange = 1L
  2379.         blue = 1L
  2380.         lt_gray = 1L
  2381.         med_green = 0L
  2382.         brown = 0L
  2383.         olive = 0L
  2384.         purple = 0L
  2385.         dk_gray = 0L
  2386.         black = 0L
  2387.     endelse
  2388.  
  2389.     colornames = [white, yellow, lavender, aqua, pink, $
  2390.         green, red, orange, blue, lt_gray, med_green, brown, $
  2391.         olive, purple, dk_gray, black]
  2392.  
  2393.     ;  Number of color indices to shift down.
  2394.     ;
  2395.     bias = FIX( (max_color ) * 0.40)
  2396.     nshift = bias 
  2397.     loadct, 23
  2398.     TVLCT, Rori, Gori, Bori, /GET
  2399.  
  2400.     r = Rori
  2401.     g = Gori
  2402.     b = Bori
  2403.  
  2404.     ncolo = !D.TABLE_SIZE    ; number of available colors
  2405.     r(0:ncolo-nshift-1) = Rori(nshift:ncolo-1)
  2406.     g(0:ncolo-nshift-1) = Gori(nshift:ncolo-1)
  2407.     b(0:ncolo-nshift-1) = Bori(nshift:ncolo-1)
  2408.  
  2409.     for i = 1, nshift do begin
  2410.         r(ncolo-nshift-1+i) = Rori(ncolo-1)
  2411.         g(ncolo-nshift-1+i) = Gori(ncolo-1)
  2412.         b(ncolo-nshift-1+i) = Bori(ncolo-1)
  2413.     endfor
  2414.  
  2415.     ;  Set up the color table.
  2416.     ;
  2417.     TVLCT, r, g, b, 0
  2418.  
  2419.     ;  Assign 16 selected colors at the top of the color table.
  2420.     ;
  2421.     opt_Init16Colors, colornames
  2422.  
  2423.     ;  Get the tips.
  2424.     ;
  2425.     sText = getTips(filepath('optimize.tip', $
  2426.         SUBDIR=['examples','demo', 'demotext']) )
  2427.  
  2428.     ;  This is the character scaling factor for
  2429.     ;  cross-platform compatibility.
  2430.     ;
  2431.     charScale = 8.0/!D.X_CH_SIZE
  2432.  
  2433.     ;  Create the widget heirarchy.
  2434.     ;
  2435.     if (ngroup EQ 0) then begin
  2436.         tlb = WIDGET_BASE(TITLE='Optimization', /COLUMN, $
  2437.             MBAR=bar_base,  $
  2438.             MAP=0, $
  2439.             /TLB_KILL_REQUEST_EVENTS, $
  2440.             TLB_FRAME_ATTR=1)
  2441.     endif else begin
  2442.         tlb = WIDGET_BASE(TITLE='Optimization', /COLUMN, $
  2443.             GROUP_LEADER=group, $
  2444.             MBAR=bar_base,  $
  2445.             MAP=0, $
  2446.             /TLB_KILL_REQUEST_EVENTS, $
  2447.             TLB_FRAME_ATTR=1)
  2448.     endelse
  2449.  
  2450.         ; Create the File|Quit buttons.
  2451.         ;
  2452.         file_menu = WIDGET_BUTTON(bar_base, VALUE='File', /MENU)
  2453.  
  2454.             quit_button = WIDGET_BUTTON(file_menu, VALUE='Quit')
  2455.  
  2456.         ; Create the Help|About buttons.
  2457.         ;
  2458.         help_menu = WIDGET_BUTTON(bar_base, VALUE='About', /HELP, /MENU)
  2459.  
  2460.             about_button = WIDGET_BUTTON(help_menu, $
  2461.                 VALUE='About Optimization')
  2462.  
  2463.  
  2464.         ;  Create the 1st base child of tlb.
  2465.         ;
  2466.         first_base = WIDGET_BASE(tlb, Column = 2)
  2467.  
  2468.             ;  Create sub base of  second base.
  2469.             ;
  2470.             subbase1 = WIDGET_BASE(first_base, Column=1)
  2471.  
  2472.                 ;  Create the selection base (surface and start).
  2473.                 ;
  2474.                 select_base = WIDGET_BASE(subbase1, /COLUMN, $
  2475.                     /BASE_ALIGN_CENTER)
  2476.  
  2477.                     ;  Create the surface pull-down menu
  2478.                     ;
  2479.                     surfaceButtonID = WIDGET_BUTTON(select_base, $
  2480.                         VALUE=' Select a Surface', Menu = 1)
  2481.  
  2482.                         surface1ID = Widget_Button(surfaceButtonID, $
  2483.                             VALUE='Surface 1')
  2484.  
  2485.                         surface2ID = Widget_Button(surfaceButtonID, $
  2486.                             VALUE='Surface 2')
  2487.  
  2488.                         surface3ID = Widget_Button(surfaceButtonID, $
  2489.                             VALUE='Surface 3')
  2490.  
  2491.                     ;  Create the staritng points pull-down menu.
  2492.                     ;
  2493.                     startButtonID = WIDGET_BUTTON(select_base, $
  2494.                         VALUE=' Select a Starting Point', Menu = 1)
  2495.  
  2496.                         point1ID = WIDGET_BUTTON(startButtonID, $
  2497.                             VALUE='Point 1')
  2498.  
  2499.                         point2ID = WIDGET_BUTTON(startButtonID, $
  2500.                             VALUE='Point 2')
  2501.  
  2502.                         point3ID = WIDGET_BUTTON(startButtonID, $
  2503.                             VALUE='Point 3')
  2504.  
  2505.                         point4ID = WIDGET_BUTTON(startButtonID, $
  2506.                             VALUE='Point 4')
  2507.  
  2508.                     ;  Create the button that activates the selections.
  2509.                     ;
  2510.                     find_button = WIDGET_BUTTON(select_base, $
  2511.                         VALUE='Find Minimum')
  2512.  
  2513.                     wResetbutton = WIDGET_BUTTON(select_base, $
  2514.                         VALUE='Reset Surface')
  2515.  
  2516.                 ;  Create a 2nd base child of subbase1
  2517.                 ;  that shows the coordinate of
  2518.                 ;  the selected starting point.
  2519.                 ;
  2520.                 show_Start_base = WIDGET_BASE(subbase1,  $
  2521.                     /COLUMN, /ALIGN_CENTER, YPAD=10)
  2522.  
  2523.                     show_start_label = WIDGET_LABEL(show_start_base, $
  2524.                         VALUE='Starting Point')
  2525.  
  2526.                     show_start_xy = WIDGET_LABEL(show_start_base, $
  2527.                         VALUE='X = -1.2    Y = 2.0')
  2528.  
  2529.                 ;  Create a 3rd base child of subbase1
  2530.                 ;  that shows the DFP results.
  2531.                 ;
  2532.                 show_RES_base = WIDGET_BASE(subbase1, /COLUMN)
  2533.  
  2534.  
  2535.                     show_RES_label = WIDGET_LABEL(show_RES_base, $
  2536.                         /ALIGN_LEFT, $
  2537.                         VALUE='Results   DFP       Powell   ')
  2538.  
  2539.                     show_RES_x = WIDGET_LABEL(show_RES_base, $
  2540.                         /ALIGN_LEFT, $
  2541.                         VALUE='X =      -1.000    -1.000    ')
  2542.  
  2543.                     show_RES_y = WIDGET_LABEL(show_RES_base, $
  2544.                         /ALIGN_LEFT, $
  2545.                         VALUE='Y =      -1.000    -1.000    ')
  2546.  
  2547.                     show_RES_z = WIDGET_LABEL(show_RES_base, $
  2548.                         /ALIGN_LEFT, $
  2549.                         VALUE='Z =      -0.000    -0.000    ')
  2550.  
  2551.                     show_RES_i = WIDGET_LABEL(show_RES_base, $
  2552.                         /ALIGN_LEFT, $
  2553.                         VALUE='#ITER =   24          26    ')
  2554.  
  2555.  
  2556.             ;   Create sub base of the first base.
  2557.             ;
  2558.             subbase2 = WIDGET_BASE(first_base)
  2559.  
  2560.                 ;  Create the drawing area.
  2561.                 ;
  2562.                 drawxsize = 0.6*screenSize(0)
  2563.                 drawysize = 0.8*drawxsize
  2564.                 drawID = Widget_DRAW(subbase2, RETAIN=2, $
  2565.                     XSIZE=drawxsize, YSIZE=drawysize)
  2566.  
  2567.         ;  Create tips texts.
  2568.         ;
  2569.         wStatusBase = WIDGET_BASE(tlb, MAP=0, /ROW)
  2570.  
  2571.             nWidgets = 2
  2572.             wText = LONARR(nWidgets)
  2573.             widTips, wStatusBase, sText.text, XSIZE=36, $
  2574.                 YSIZE=3, NWIDGETS=nWidgets, wText
  2575.  
  2576.     ;  All the widgets has been created.
  2577.  
  2578.     ;  Realize the top-level base.
  2579.     ;
  2580.     WIDGET_CONTROL, tlb, /REALIZE
  2581.  
  2582.     ;  Returns the top level base to the APPTLB keyword.
  2583.     ;
  2584.     appTLB = tlb
  2585.  
  2586.     ;  Size the tips widgets.
  2587.     ;
  2588.     sizeTips, tlb, wText, wStatusBase
  2589.  
  2590.     ; Get the drawing window index
  2591.     ;
  2592.     WIDGET_CONTROL, drawID, GET_VALUE=drawWindowID
  2593.  
  2594.     ;  Create the 3 surfaces.
  2595.     ;  Surface index 1 first (banana function
  2596.     ;
  2597.     surfaceIndex = 1
  2598.     dim_array = 33   ; dimension of the array
  2599.     zdata1 = DBLARR(dim_array, dim_array)
  2600.  
  2601.     xdata1 = FINDGEN(dim_array)/8.0 -2.0
  2602.     ydata1 = FINDGEN(dim_array)/8.0 -1.0
  2603.  
  2604.     ; Create the 2-D arrays of the mesh points
  2605.     ;
  2606.     for i = 0, dim_array-1 do begin
  2607.         for j = 0 , dim_array-1 do begin
  2608.             x = [xdata1(i), ydata1(j)]
  2609.             zdata1(i,j) = FUNCOP(x, surfaceIndex)
  2610.         endfor
  2611.     endfor
  2612.  
  2613.     ;  Create the surface 2 and surface 3 data SETS here.
  2614.     ;
  2615.     surfaceindex = 2
  2616.     dim_array = 33
  2617.     zdata2 = DBLARR(dim_array, dim_array)
  2618.     xdata2 = FINDGEN(dim_array)/8.0 -2.0
  2619.     ydata2 = FINDGEN(dim_array)/8.0 -2.0
  2620.  
  2621.     ;  Create the 2-D arrays of the mesh points.
  2622.     ;
  2623.     for i = 0, dim_array-1 do begin
  2624.         for j = 0 , dim_array-1 do begin
  2625.             x = [xdata2(i), ydata2(j)]
  2626.             zdata2(i,j) = FUNCOP(x, surfaceIndex)
  2627.         endfor
  2628.     endfor
  2629.  
  2630.     ; Create the surface 3  data SETS here
  2631.     ;
  2632.     surfaceindex = 3
  2633.  
  2634.     dim_array = 17
  2635.     zdata3 = DBLARR(dim_array, dim_array)
  2636.  
  2637.     xdata3 = FINDGEN(dim_array)/8.0 -1.0
  2638.     ydata3 = FINDGEN(dim_array)/8.0 -1.0
  2639.  
  2640.     ; Create the 2-D arrays of the mesh points
  2641.     ;
  2642.     for i = 0, dim_array-1 do begin
  2643.         for j = 0 , dim_array-1 do begin
  2644.             x = [xdata3(i), ydata3(j)]
  2645.             zdata3(i,j) = FUNCOP(x, surfaceIndex)
  2646.         endfor
  2647.     endfor
  2648.  
  2649.     ;  Draw the default surface (surface 1) on a newly
  2650.     ;  created pixmap.
  2651.     ;
  2652.     WINDOW, /FREE, XSIZE=drawxsize, YSIZE=drawysize, /PIXMAP
  2653.     pixmap1ID = !d.Window
  2654.  
  2655.     WSET, pixmap1ID
  2656.  
  2657.     ;  Plot the surface 1 on the pixmap
  2658.     ;  and save the 3d tranformation
  2659.     ;
  2660.     if (max_color EQ 16) then begin
  2661.         topc = 15
  2662.     endif else if (max_color GE 17) then begin
  2663.         topc = max_color - 16
  2664.     endif else begin
  2665.         topc = 1
  2666.     endelse
  2667.  
  2668.     ;  Draw the shaded surface of the Banana function without axes.
  2669.     ;
  2670.     SHADE_SURF, zdata1, xdata1, ydata1, /SAVE, $
  2671.         SHADES= BYTSCL(zdata1, TOP=topc), $
  2672.         BACKGROUND=black, AZ=10, AX= 60, $
  2673.         XSTYLE=4, YSTYLE=4, ZSTYLE=4, $
  2674.         POSITION=[0.1, 0.20, 0.95, 1.]
  2675.  
  2676.     ; Draw the mesh surface of the Banana function without axes
  2677.     ;
  2678.     SURFACE, zdata1, xdata1, ydata1, /SAVE,/NOERASE, $
  2679.         COLOR = 31,AZ=10, AX = 60, $
  2680.         XSTYLE=4, YSTYLE=4, ZSTYLE=4, $
  2681.         BACKGROUND=black, $
  2682.         POSITION=[0.1, 0.20, 0.95, 1.]
  2683.  
  2684.     saved3d1 = !p.t
  2685.     !P.T = 0
  2686.  
  2687.     ; Now create the pixmap for the 2nd surface
  2688.     ;
  2689.  
  2690.     ; Save the current color table
  2691.     ;
  2692.     TVLCT, savedR2, savedG2, savedB2, /GET
  2693.     colorTable2=[[savedR2], [savedG2], [savedB2]]
  2694.  
  2695.     WINDOW, /FREE, XSIZE=drawxsize, YSIZE=drawysize, /PIXMAP
  2696.     pixmap2ID = !d.Window
  2697.  
  2698.     !p.T = 0
  2699.     WSET, pixmap2ID
  2700.     TVLCT, Rori, Gori, Bori,0
  2701.  
  2702.     SHADE_SURF, zdata2, xdata2, ydata2, /SAVE, $
  2703.         SHADES=BYTSCL(zdata2, TOP=topc-nshift), $
  2704.         BACKGROUND=black, AZ=30, AX=30, $
  2705.         XRANGE=[-2,2], $
  2706.         YRANGE=[-2,2], $
  2707.         ZRANGE=[-7,7],$
  2708.         XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  2709.         POSITION=[0.05, 0.05, 0.95, 1.]
  2710.  
  2711.     SURFACE, zdata2, xdata2, ydata2, /SAVE,/NOERASE, $
  2712.         BACKGROUND=black, AZ=30, AX=30, $
  2713.         XRANGE=[-2,2], $
  2714.         YRANGE=[-2,2], $
  2715.         ZRANGE=[-7,7],$
  2716.         COLOR=31, $
  2717.         XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  2718.         POSITION=[0.05, 0.05, 0.95, 1.]
  2719.  
  2720.     saved3d2 = !p.t
  2721.     TVLCT,colortable2
  2722.  
  2723.     ;  Now create the pixmap for the 3rd surface.
  2724.     ;
  2725.     ;  Save the current color table.
  2726.     ;
  2727.     TVLCT, savedR3, savedG3, savedB3,/GET
  2728.     colorTable3=[[savedR3], [savedG3], [savedB3]]
  2729.  
  2730.     WINDOW, /FREE, XSIZE=drawxsize, YSIZE=drawysize, /PIXMAP
  2731.     pixmap3ID = !d.Window
  2732.  
  2733.     !p.T = 0
  2734.     WSET, pixmap3ID
  2735.     TVLCT, Rori, Gori, Bori,0
  2736.     SHADE_SURF, zdata3, xdata3, ydata3, /SAVE, $
  2737.         shades = BYTSCL(zdata3, TOP = topc-nshift), $
  2738.         BACKGROUND=black, AZ=30, AX=30, $
  2739.         XRANGE=[-1,1], $
  2740.         YRANGE=[-1,1], $
  2741.         ZRANGE=[-1,1],$
  2742.         XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  2743.         POSITION=[0.05, 0.05, 0.95, 1.]
  2744.  
  2745.     SURFACE, zdata3, xdata3, ydata3, /SAVE,/NOERASE, $
  2746.         BACKGROUND=black, AZ=30, AX=30, $
  2747.         XRANGE=[-1,1], $
  2748.         YRANGE=[-1,1], $
  2749.         ZRANGE=[-1,1],$
  2750.         COLOR=31, $
  2751.         XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  2752.         POSITION=[0.05, 0.05, 0.95, 1.]
  2753.  
  2754.     saved3d3 = !p.t
  2755.  
  2756.     WSET, drawWindowID
  2757.     !P.T = 0
  2758.     TVLCT,colortable3
  2759.  
  2760.     ;  Initialize the  starting points.
  2761.     ;
  2762.     ;  1st index : surface index (1,2,3)
  2763.     ;  2nd index : point index (1,2,3,4)
  2764.     ;  3rd index ; x, y, z data coordinates
  2765.     ;
  2766.     startPoints = DBLARR(4,5,4)
  2767.     startPointsdev = DBLARR(4,5,4)
  2768.     opt_InitPoint, startpoints
  2769.  
  2770.     ;  Convert the surface 1 data into device coordinates.
  2771.     ;
  2772.     !P.T = 0
  2773.     !P.T = saved3d1
  2774.  
  2775.     d = DBLARR(3,1)
  2776.  
  2777.     ;  Here we create an empty surface in order to
  2778.     ;  establish the data system coordinates.
  2779.     ;  Then, we are able to convert the data into device coordinates.
  2780.     ;
  2781.     SURFACE, zdata1, xdata1, ydata1,/nodata,  $
  2782.         COLOR = 31,AZ=10, AX = 60, $
  2783.         XSTYLE=4, YSTYLE=4, ZSTYLE=4, $
  2784.         BACKGROUND=black, $
  2785.         POSITION=[0.1, 0.20, 0.95, 1.]
  2786.  
  2787.     for i = 1, 4 do begin
  2788.         d = CONVERT_COORD(startPoints(1, i, 1), $
  2789.             startPoints(1, i, 2), $
  2790.             startPoints(1, i, 3), $
  2791.             /TO_DEVICE, /T3D)
  2792.         startPointsdev(1, i, 1:3) = d(0:2, 0)
  2793.     endfor
  2794.  
  2795.     ; Convert the surface 2 data into device coordinates
  2796.     ;
  2797.     !P.T = 0
  2798.     !P.T = saved3d2
  2799.  
  2800.     SURFACE, zdata2, xdata2, $
  2801.         ydata2,  /NODATA, $
  2802.         BACKGROUND=black, AZ=30, AX=30,$
  2803.         XRANGE=[-2,2], $
  2804.         YRANGE=[-2,2], $
  2805.         ZRANGE=[-7,7],$
  2806.         XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  2807.         POSITION=[0.05, 0.05, 0.95, 1.]
  2808.  
  2809.     d = DBLARR(3,1)
  2810.     for i = 1, 4 do begin
  2811.         d = CONVERT_COORD(startPoints(2, i, 1), $
  2812.             startPoints(2, i, 2), $
  2813.             startPoints(2, i, 3), $
  2814.             /TO_DEVICE, /T3D)
  2815.         startPointsdev(2, i, 1:3) = d(0:2, 0)
  2816.     endfor
  2817.  
  2818.     ; Convert the surface 3 data into device coordinates
  2819.     ;
  2820.     !P.T = 0
  2821.     !P.T = saved3d3
  2822.  
  2823.     SURFACE, zdata3, xdata3, ydata3, /SAVE,/NODATA, $
  2824.         BACKGROUND=black, AZ=30, AX=30,$
  2825.         XRANGE=[-1,1], $
  2826.         YRANGE=[-1,1], $
  2827.         ZRANGE=[-1,1],$
  2828.         XSTYLE=6, YSTYLE=6, ZSTYLE=6, $
  2829.         POSITION=[0.05, 0.05, 0.95, 1.]
  2830.  
  2831.     d = DBLARR(3,1)
  2832.     for i = 1, 4 do begin
  2833.         d = CONVERT_COORD(startPoints(3, i, 1), $
  2834.             startPoints(3, i, 2), $
  2835.             startPoints(3, i, 3), $
  2836.             /TO_DEVICE, /T3D)
  2837.         startPointsdev(3, i, 1:3) = d(0:2, 0)
  2838.     endfor
  2839.  
  2840.     !P.T = 0
  2841.  
  2842.     ;  Point Flag  0 =  No starting point has been selected
  2843.     ;              1 =  A starting point has been selected
  2844.     ;
  2845.     pointFlag = 0
  2846.  
  2847.     ;  Desensitize the Find minimum button.
  2848.     ;
  2849.     WIDGET_CONTROL, find_button, sensitive = 0
  2850.     WIDGET_CONTROL, surface1ID, sensitive = 0
  2851.  
  2852.     ;  Define the info structure.
  2853.     ;
  2854.     info = { $
  2855.         SurfaceButtonID: surfacebuttonID, $   ; Surface button
  2856.         Surface1ID: surface1ID, $             ; Surfaces 1,2,3 buttons
  2857.         Surface2ID: surface2ID, $
  2858.         Surface3ID: surface3ID, $
  2859.         WResetButton: wResetButton, $         ; Reset surface button
  2860.         SurfaceIndex: 1, $                    ; Default surface is no.1
  2861.         StartButtonID: startbuttonID, $       ; Starting Point buttons
  2862.         AboutButtonID: about_button, $        ; Information on the applet
  2863.         Point1ID: point1ID, $
  2864.         Point2ID: point2ID, $
  2865.         Point3ID: point3ID, $
  2866.         Point4ID: point4ID, $
  2867.         FindButton: find_button, $            ; Find button
  2868.         StartLbL: show_start_xy, $            ; Starting point label
  2869.         ResultLbLx : show_res_x, $            ; Results labels
  2870.         ResultLbLy : show_res_y, $
  2871.         ResultLbLz : show_res_z, $
  2872.         ResultLbLi : show_res_i, $            ; Iterations
  2873.         Quitbutton: quit_button, $            ; Quit button in file menu
  2874.         Pixmap1ID: pixmap1ID, $               ; Pixmaps ID's
  2875.         Pixmap2ID: pixmap2ID, $
  2876.         Pixmap3ID :pixmap3ID, $
  2877.         DrawWinID: drawwindowID, $            ; Drawing window
  2878.         Xdata1: xdata1, $                     ; Surface data sets
  2879.         Ydata1: ydata1, $
  2880.         Zdata1: zdata1, $
  2881.         Xdata2: xdata2, $
  2882.         Ydata2: ydata2, $
  2883.         Zdata2: zdata2, $
  2884.         Xdata3: xdata3, $
  2885.         Ydata3: ydata3, $
  2886.         Zdata3: zdata3, $
  2887.         ColorTable: colortable, $             ; Color table to restore
  2888.         PointIndex: 0, $                      ; Selected point
  2889.         PointFlag: pointFlag, $               ; If point selected(0 = no)
  2890.         StartPoints: startpoints, $           ; Starting points in data coor.
  2891.         StartPointsdev: startpointsdev, $     ; Start. points in device coor.
  2892.         Saved3d1: saved3d1, $                 ; 3-D transformation surface 1
  2893.         Saved3d2: saved3d2, $                 ; 3-D transformation surface 2
  2894.         Saved3d3: saved3d3, $                 ; 3-D transformation surface 3
  2895.         drawxsize: drawxsize, $               ; Drawing window size
  2896.         drawysize: drawysize, $
  2897.         Red: red, $                           ; Colors indices
  2898.         Yellow: yellow, $
  2899.         Lavender: lavender, $
  2900.         Aqua: aqua, $
  2901.         Pink: pink, $
  2902.         White: white, $
  2903.         Green: green, $
  2904.         Orange: orange, $
  2905.         Blue: blue, $
  2906.         Lt_gray: lt_gray, $
  2907.         Med_green: med_green, $
  2908.         Brown: brown, $
  2909.         Olive: olive, $
  2910.         Purple: purple, $
  2911.         Dk_gray: dk_gray, $
  2912.         Black: black, $
  2913.         FontFlag : fontflag, $
  2914.         CharScale: charScale, $                ; Character scaling factor
  2915.         groupBase: groupBase $                 ; Base of Group Leader
  2916.     }
  2917.  
  2918.     ;  Reset the result label to default.
  2919.     ;
  2920.     opt_ResetResult, info.ResultLBLx, info.ResultLBLy, $
  2921.         info.ResultLBLz, info.ResultLBLi
  2922.  
  2923.     ;  Set the info structure.
  2924.     ;
  2925.     WIDGET_CONTROL, tlb, SET_UVALUE=info, /NO_COPY
  2926.  
  2927.  
  2928.     ;  Copy surface 1 into the drawing window.
  2929.     ;
  2930.     WSET, drawWindowID
  2931.     SURFACE, zdata1, xdata1, ydata1, /NODATA,  $
  2932.         COLOR=31, AZ=10, AX=60, $
  2933.         XStyle=4, YStyle=4, ZStyle=4, $
  2934.         POSITION=[0.1, 0.20, 0.95, 1.],background = black
  2935.  
  2936.     Device, Copy=[0, 0, drawxsize, drawysize, 0, 0, pixmap1ID]
  2937.  
  2938.     ;  Plot the 4 starting points into the drawing window.
  2939.     ;
  2940.     surfaceIndex = 1
  2941.     pointnumber = 1
  2942.     saved3dt = saved3d1
  2943.     opt_PlotPoint, surfaceIndex, pointnumber, charscale,red, $
  2944.         saved3dt, startPoints, startPointsdev, FontFlag,  /all
  2945.  
  2946.  
  2947.     ;  Destroy the starting up window.
  2948.     ;
  2949.     WIDGET_CONTROL, drawbase, /DESTROY
  2950.  
  2951.     ;  Map the top level base.
  2952.     ;
  2953.     WIDGET_CONTROL, tlb, MAP=1
  2954.  
  2955.     ;  Register with the XManager.
  2956.     ;
  2957.     XMANAGER, "D_OPTIMIZE", tlb, Event_Handler="OPTIMIZE_EVENT",$
  2958.         CLEANUP='OPTIMIZE_CLEANUP', $
  2959.         /NO_BLOCK
  2960.  
  2961. end   ; OF D_OPTIMIZE.PRO
  2962.