function SCAT_LOG, X, P

                 RETURN, P[0] *  (  P[1]  -  ALOG10 (X-1) )
END

;---  pro memu

PRO MENU, BUT_DESC, OPTION, TITLE=TITLE, XOFFSET=XOFFSET, YOFFSET=YOFFSET, $
      ROW=ROW, BUTSIZE=BUTSIZE, BASE=BASE, BUTTONS=BUTTONS, FONT=FONT, $
      frame=frame, ypad=ypad, no_release = no_release, $
      scroll=scroll, space=space, xpad=xpad, exclusive=exclusive, $
      uvalue=uvalue, x_scroll_size=x_scroll_size, $
      y_scroll_size=y_scroll_size


parent = 0
if (not keyword_set(TITLE)) then TITLE = 'MENU'
s = size(but_desc)
value_type = s(s(0) + 1)
if ((value_type ne 1) and (value_type ne 7)) then $
    message, 'values must be a string vector or 3-d byte array.`
  if (value_type eq 1) then begin
    if (s(0) ne 3) then message, 'type byte values must be 3-d'
    n_buttons = s(3)
  endif else begin
    n_buttons = n_elements(but_desc)
endelse

; sort out the keywords
if (not keyword_set(butsize)) then butsize=100
if (not keyword_set(xoffset)) then xoffset=5
if (not keyword_set(yoffset)) then yoffset=25
if (not keyword_set(row)) then begin
    row=0
    column=1
  endif else begin
    row=1
    column=0
endelse
if (not keyword_set(font)) then font = ''
if (not keyword_set(exclusive)) then exclusive=0
if (not keyword_set(nonexclusive)) then nonexclusive=0
if (keyword_set(scroll) or keyword_set(x_scroll_size) or $
      keyword_set(y_scroll_size)) then begin
     scroll = 1;
     if (not keyword_set(x_scroll_size)) then x_scroll_size=0
     if (not keyword_set(y_scroll_size)) then y_scroll_size=0
     endif else begin
     scroll=0
endelse
if (not keyword_set(frame)) then frame = 0
if (not keyword_set(space)) then space = 0
if (not keyword_set(xpad)) then xpad = 0
if (not keyword_set(ypad)) then ypad = 0
if (not keyword_set(uvalue)) then begin
    uvalue=lindgen(n_buttons)
  endif else begin
    s = size(uvalue)
    if (s(s(0) + 2) ne n_buttons) then $
      message, 'uvalue must have the same number of elements as values'
endelse

; create the base
  if (parent eq 0) then begin
    if (scroll) then $
      base = widget_base(column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, scroll=scroll, $
      space=space, xpad=xpad, ypad=ypad, x_scroll_size=x_scroll_size, $
      y_scroll_size=y_scroll_size, title = title) $
    else $
      base = widget_base(column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, $
      space=space, xpad=xpad, ypad=ypad, title = title)
  endif else begin
    if (keyword_set(title)) then begin
      theparent = widget_base(parent, /column, /frame)
      thelabel = widget_label(theparent, value = title)
    endif else theparent = parent
    if (scroll) then $
      base = widget_base(theparent, column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, scroll=scroll, $
      space=space, xpad=xpad, ypad=ypad, x_scroll_size=x_scroll_size, $
      y_scroll_size=y_scroll_size) $
    else $
      base = widget_base(theparent, column=column, exclusive=exclusive, $
      frame=frame, nonexclusive=nonexclusive, row=row, $
      space=space, xpad=xpad, ypad=ypad)
  endelse

; create the buttons
buttons = lindgen(n_buttons)
if (value_type eq 1) then begin
    for i = 0, n_buttons-1 do $
      buttons(i) = widget_button(base, value=but_desc(*, *, i), $
       no_release = no_release, uvalue=uvalue(i), $
       font = font, xsize=butsize)
endif else begin
    for i = 0, n_buttons-1 do $
      buttons(i) = widget_button(base, value=but_desc(i), $
       no_release = no_release, uvalue=uvalue(i), $
       font = font, xsize=butsize)
endelse

; realizacja na ekranie

widget_control, /realize, base, tlb_set_xoffset=xoffset, tlb_set_yoffset=yoffset
event=widget_event(base)
kl= fix ( where (buttons eq event.id) + 1 )
widget_control, base, /destroy
option=kl(0)
end



;---  end of pro menu   ----------------------------------------

;+
; NAME:
;   MPFITFUN
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE:
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Perform Levenberg-Marquardt least-squares fit to IDL function
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   parms = MPFITFUN(MYFUNCT, X, Y, ERR, start_params, ...)
;
; DESCRIPTION:
;
;  MPFITFUN fits a user-supplied model -- in the form of an IDL
;  function -- to a set of user-supplied data.  MPFITFUN calls
;  MPFIT, the MINPACK-1 least-squares minimizer, to do the main
;  work.
;
;  Given the data and their uncertainties, MPFITFUN finds the best set
;  of model parameters which match the data (in a least-squares
;  sense) and returns them in an array.
;
;  The user must supply the following items:
;   - An array of independent variable values ("X").
;   - An array of "measured" *dependent* variable values ("Y").
;   - An array of "measured" 1-sigma uncertainty values ("ERR").
;   - The name of an IDL function which computes Y given X ("MYFUNCT").
;   - Starting guesses for all of the parameters ("START_PARAMS").
;
;  There are very few restrictions placed on X, Y or MYFUNCT.  Simply
;  put, MYFUNCT must map the "X" values into "Y" values given the
;  model parameters.  The "X" values may represent any independent
;  variable (not just Cartesian X), and indeed may be multidimensional
;  themselves.  For example, in the application of image fitting, X
;  may be a 2xN array of image positions.
;
;  Data values of NaN or Infinity for "Y", "ERR" or "WEIGHTS" will be
;  ignored as missing data if the NAN keyword is set.  Otherwise, they
;  may cause the fitting loop to halt with an error message.  Note
;  that the fit will still halt if the model function, or its
;  derivatives, produces infinite or NaN values.
;
;  MPFITFUN carefully avoids passing large arrays where possible to
;  improve performance.
;
;  See below for an example of usage.
;
; USER FUNCTION
;
;  The user must define a function which returns the model value.  For
;  applications which use finite-difference derivatives -- the default
;  -- the user function should be declared in the following way:
;
;    FUNCTION MYFUNCT, X, P
;     ; The independent variable is X
;     ; Parameter values are passed in "P"
;     YMOD = ... computed model values at X ...
;     return, YMOD
;    END
;
;  The returned array YMOD must have the same dimensions and type as
;  the "measured" Y values.
;
;  User functions may also indicate a fatal error condition
;  using the ERROR_CODE common block variable, as described
;  below under the MPFIT_ERROR common block definition.
;
;  See the discussion under "EXPLICIT DERIVATIVES" and AUTODERIVATIVE
;  in MPFIT.PRO if you wish to compute the derivatives for yourself.
;  AUTODERIVATIVE is accepted by MPFITFUN and passed directly to
;  MPFIT.  The user function must accept one additional parameter, DP,
;  which contains the derivative of the user function with respect to
;  each parameter at each data point, as described in MPFIT.PRO.
;
; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD
;
;  The behavior of MPFIT can be modified with respect to each
;  parameter to be fitted.  A parameter value can be fixed; simple
;  boundary constraints can be imposed; limitations on the parameter
;  changes can be imposed; properties of the automatic derivative can
;  be modified; and parameters can be tied to one another.
;
;  These properties are governed by the PARINFO structure, which is
;  passed as a keyword parameter to MPFIT.
;
;  PARINFO should be an array of structures, one for each parameter.
;  Each parameter is associated with one element of the array, in
;  numerical order.  The structure can have the following entries
;  (none are required):
;
;     .VALUE - the starting parameter value (but see the START_PARAMS
;              parameter for more information).
;
;     .FIXED - a boolean value, whether the parameter is to be held
;              fixed or not.  Fixed parameters are not varied by
;              MPFIT, but are passed on to MYFUNCT for evaluation.
;
;     .LIMITED - a two-element boolean array.  If the first/second
;                element is set, then the parameter is bounded on the
;                lower/upper side.  A parameter can be bounded on both
;                sides.  Both LIMITED and LIMITS must be given
;                together.
;
;     .LIMITS - a two-element float or double array.  Gives the
;               parameter limits on the lower and upper sides,
;               respectively.  Zero, one or two of these values can be
;               set, depending on the values of LIMITED.  Both LIMITED
;               and LIMITS must be given together.
;
;     .PARNAME - a string, giving the name of the parameter.  The
;                fitting code of MPFIT does not use this tag in any
;                way.  However, the default ITERPROC will print the
;                parameter name if available.
;
;     .STEP - the step size to be used in calculating the numerical
;             derivatives.  If set to zero, then the step size is
;             computed automatically.  Ignored when AUTODERIVATIVE=0.
;             This value is superceded by the RELSTEP value.
;
;     .RELSTEP - the *relative* step size to be used in calculating
;                the numerical derivatives.  This number is the
;                fractional size of the step, compared to the
;                parameter value.  This value supercedes the STEP
;                setting.  If the parameter is zero, then a default
;                step size is chosen.
;
;     .MPSIDE - the sidedness of the finite difference when computing
;               numerical derivatives.  This field can take four
;               values:
;
;                  0 - one-sided derivative computed automatically
;                  1 - one-sided derivative (f(x+h) - f(x)  )/h
;                 -1 - one-sided derivative (f(x)   - f(x-h))/h
;                  2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)
;
;              Where H is the STEP parameter described above.  The
;              "automatic" one-sided derivative method will chose a
;              direction for the finite difference which does not
;              violate any constraints.  The other methods do not
;              perform this check.  The two-sided method is in
;              principle more precise, but requires twice as many
;              function evaluations.  Default: 0.
;
;     .MPMAXSTEP - the maximum change to be made in the parameter
;                  value.  During the fitting process, the parameter
;                  will never be changed by more than this value in
;                  one iteration.
;
;                  A value of 0 indicates no maximum.  Default: 0.
;
;     .TIED - a string expression which "ties" the parameter to other
;             free or fixed parameters as an equality constraint.  Any
;             expression involving constants and the parameter array P
;             are permitted.
;             Example: if parameter 2 is always to be twice parameter
;             1 then use the following: parinfo[2].tied = '2 * P[1]'.
;             Since they are totally constrained, tied parameters are
;             considered to be fixed; no errors are computed for them.
;             [ NOTE: the PARNAME can't be used in a TIED expression. ]
;
;     .MPPRINT - if set to 1, then the default ITERPROC will print the
;                parameter value.  If set to 0, the parameter value
;                will not be printed.  This tag can be used to
;                selectively print only a few parameter values out of
;                many.  Default: 1 (all parameters printed)
;
;     .MPFORMAT - IDL format string to print the parameter within
;                 ITERPROC.  Default: '(G20.6)'  (An empty string will
;                 also use the default.)
;
;  Future modifications to the PARINFO structure, if any, will involve
;  adding structure tags beginning with the two letters "MP".
;  Therefore programmers are urged to avoid using tags starting with
;  "MP", but otherwise they are free to include their own fields
;  within the PARINFO structure, which will be ignored by MPFIT.
;
;  PARINFO Example:
;  parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $
;                       limits:[0.D,0]}, 5)
;  parinfo[0].fixed = 1
;  parinfo[4].limited[0] = 1
;  parinfo[4].limits[0]  = 50.D
;  parinfo[*].value = [5.7D, 2.2, 500., 1.5, 2000.]
;
;  A total of 5 parameters, with starting values of 5.7,
;  2.2, 500, 1.5, and 2000 are given.  The first parameter
;  is fixed at a value of 5.7, and the last parameter is
;  constrained to be above 50.
;
; COMPATIBILITY
;
;  This function is designed to work with IDL 5.0 or greater.
;
;  Because TIED parameters rely on the EXECUTE() function, they cannot
;  be used with the free version of the IDL Virtual Machine.
;
;
; INPUTS:
;   MYFUNCT - a string variable containing the name of an IDL function.
;             This function computes the "model" Y values given the
;             X values and model parameters, as desribed above.
;
;   X - Array of independent variable values.
;
;   Y - Array of "measured" dependent variable values.  Y should have
;       the same data type as X.  The function MYFUNCT should map
;       X->Y.
;         NOTE: the following special cases apply:
;                * if Y is NaN or Infinite, and the NAN keyword is
;                  set, then the corresponding data point is ignored
;
;   ERR - Array of "measured" 1-sigma uncertainties.  ERR should have
;         the same data type as Y.  ERR is ignored if the WEIGHTS
;         keyword is specified.
;         NOTE: the following special cases apply:
;                * if ERR is zero, then the corresponding data point
;                  is ignored
;                * if ERR is NaN or Infinite, and the NAN keyword is
;                  set, then the corresponding data point is ignored
;                * if ERR is negative, then the absolute value of
;                  ERR is used.
;
;   START_PARAMS - An array of starting values for each of the
;                  parameters of the model.  The number of parameters
;                  should be fewer than the number of measurements.
;                  Also, the parameters should have the same data type
;                  as the measurements (double is preferred).
;
;                  This parameter is optional if the PARINFO keyword
;                  is used (see MPFIT).  The PARINFO keyword provides
;                  a mechanism to fix or constrain individual
;                  parameters.  If both START_PARAMS and PARINFO are
;                  passed, then the starting *value* is taken from
;                  START_PARAMS, but the *constraints* are taken from
;                  PARINFO.
;
;
; RETURNS:
;
;   Returns the array of best-fit parameters.
;
;
; KEYWORD PARAMETERS:
;
;   BESTNORM - the value of the summed squared residuals for the
;              returned parameter values.
;
;   COVAR - the covariance matrix for the set of parameters returned
;           by MPFIT.  The matrix is NxN where N is the number of
;           parameters.  The square root of the diagonal elements
;           gives the formal 1-sigma statistical errors on the
;           parameters IF errors were treated "properly" in MYFUNC.
;           Parameter errors are also returned in PERROR.
;
;           To compute the correlation matrix, PCOR, use this:
;           IDL> PCOR = COV * 0
;           IDL> FOR i = 0, n-1 DO FOR j = 0, n-1 DO $
;                PCOR[i,j] = COV[i,j]/sqrt(COV[i,i]*COV[j,j])
;
;           If NOCOVAR is set or MPFIT terminated abnormally, then
;           COVAR is set to a scalar with value !VALUES.D_NAN.
;
;   CASH - when set, the fit statistic is changed to a derivative of
;          the CASH statistic.  The model function must be strictly
;          positive. WARNING: this option is incomplete and untested.
;
;   DOF - number of degrees of freedom, computed as
;             DOF = N_ELEMENTS(DEVIATES) - NFREE
;         Note that this doesn't account for pegged parameters (see
;         NPEGGED).
;
;   ERRMSG - a string error or warning message is returned.
;
;   FTOL - a nonnegative input variable. Termination occurs when both
;          the actual and predicted relative reductions in the sum of
;          squares are at most FTOL (and STATUS is accordingly set to
;          1 or 3).  Therefore, FTOL measures the relative error
;          desired in the sum of squares.  Default: 1D-10
;
;   FUNCTARGS - A structure which contains the parameters to be passed
;               to the user-supplied function specified by MYFUNCT via
;               the _EXTRA mechanism.  This is the way you can pass
;               additional data to your user-supplied function without
;               using common blocks.
;
;               By default, no extra parameters are passed to the
;               user-supplied function.
;
;   GTOL - a nonnegative input variable. Termination occurs when the
;          cosine of the angle between fvec and any column of the
;          jacobian is at most GTOL in absolute value (and STATUS is
;          accordingly set to 4). Therefore, GTOL measures the
;          orthogonality desired between the function vector and the
;          columns of the jacobian.  Default: 1D-10
;
;   ITERARGS - The keyword arguments to be passed to ITERPROC via the
;              _EXTRA mechanism.  This should be a structure, and is
;              similar in operation to FUNCTARGS.
;              Default: no arguments are passed.
;
;   ITERPROC - The name of a procedure to be called upon each NPRINT
;              iteration of the MPFIT routine.  It should be declared
;              in the following way:
;
;              PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $
;                PARINFO=parinfo, QUIET=quiet, ...
;                ; perform custom iteration update
;              END
;
;              ITERPROC must either accept all three keyword
;              parameters (FUNCTARGS, PARINFO and QUIET), or at least
;              accept them via the _EXTRA keyword.
;
;              MYFUNCT is the user-supplied function to be minimized,
;              P is the current set of model parameters, ITER is the
;              iteration number, and FUNCTARGS are the arguments to be
;              passed to MYFUNCT.  FNORM should be the
;              chi-squared value.  QUIET is set when no textual output
;              should be printed.  See below for documentation of
;              PARINFO.
;
;              In implementation, ITERPROC can perform updates to the
;              terminal or graphical user interface, to provide
;              feedback while the fit proceeds.  If the fit is to be
;              stopped for any reason, then ITERPROC should set the
;              common block variable ERROR_CODE to negative value (see
;              MPFIT_ERROR common block below).  In principle,
;              ITERPROC should probably not modify the parameter
;              values, because it may interfere with the algorithm's
;              stability.  In practice it is allowed.
;
;              Default: an internal routine is used to print the
;                       parameter values.
;
;   MAXITER - The maximum number of iterations to perform.  If the
;             number is exceeded, then the STATUS value is set to 5
;             and MPFIT returns.
;             Default: 200 iterations
;
;   NAN - ignore infinite or NaN values in the Y, ERR or WEIGHTS
;         parameters.  These values will be treated as missing data.
;         However, the fit will still halt with an error condition
;         if the model function becomes infinite.
;
;   NFEV - the number of MYFUNCT function evaluations performed.
;
;   NFREE - the number of free parameters in the fit.  This includes
;           parameters which are not FIXED and not TIED, but it does
;           include parameters which are pegged at LIMITS.
;
;   NITER - the number of iterations completed.
;
;   NOCOVAR - set this keyword to prevent the calculation of the
;             covariance matrix before returning (see COVAR)
;
;   NPEGGED - the number of free parameters which are pegged at a
;             LIMIT.
;
;   NPRINT - The frequency with which ITERPROC is called.  A value of
;            1 indicates that ITERPROC is called with every iteration,
;            while 2 indicates every other iteration, etc.  Note that
;            several Levenberg-Marquardt attempts can be made in a
;            single iteration.
;            Default value: 1
;
;   PARINFO - Provides a mechanism for more sophisticated constraints
;             to be placed on parameter values.  When PARINFO is not
;             passed, then it is assumed that all parameters are free
;             and unconstrained.  Values in PARINFO are never
;             modified during a call to MPFIT.
;
;             See description above for the structure of PARINFO.
;
;             Default value:  all parameters are free and unconstrained.
;
;   PERROR - The formal 1-sigma errors in each parameter, computed
;            from the covariance matrix.  If a parameter is held
;            fixed, or if it touches a boundary, then the error is
;            reported as zero.
;
;            If the fit is unweighted (i.e. no errors were given, or
;            the weights were uniformly set to unity), then PERROR
;            will probably not represent the true parameter
;            uncertainties.
;
;            *If* you can assume that the true reduced chi-squared
;            value is unity -- meaning that the fit is implicitly
;            assumed to be of good quality -- then the estimated
;            parameter uncertainties can be computed by scaling PERROR
;            by the measured chi-squared value.
;
;              DOF     = N_ELEMENTS(X) - N_ELEMENTS(PARMS) ; deg of freedom
;              PCERROR = PERROR * SQRT(BESTNORM / DOF)   ; scaled uncertainties
;
;   QUIET - set this keyword when no textual output should be printed
;           by MPFIT
;
;   STATUS - an integer status code is returned.  All values other
;            than zero can represent success.  It can have one of the
;            following values:
;
;      0  improper input parameters.
;
;      1  both actual and predicted relative reductions
;         in the sum of squares are at most FTOL.
;
;      2  relative error between two consecutive iterates
;         is at most XTOL
;
;      3  conditions for STATUS = 1 and STATUS = 2 both hold.
;
;      4  the cosine of the angle between fvec and any
;         column of the jacobian is at most GTOL in
;         absolute value.
;
;      5  the maximum number of iterations has been reached
;
;      6  FTOL is too small. no further reduction in
;         the sum of squares is possible.
;
;      7  XTOL is too small. no further improvement in
;         the approximate solution x is possible.
;
;      8  GTOL is too small. fvec is orthogonal to the
;         columns of the jacobian to machine precision.
;
;   WEIGHTS - Array of weights to be used in calculating the
;             chi-squared value.  If WEIGHTS is specified then the ERR
;             parameter is ignored.  The chi-squared value is computed
;             as follows:
;
;                CHISQ = TOTAL( (Y-MYFUNCT(X,P))^2 * ABS(WEIGHTS) )
;
;             Here are common values of WEIGHTS for standard weightings:
;
;                1D/ERR^2 - Normal weighting (ERR is the measurement error)
;                1D/Y     - Poisson weighting (counting statistics)
;                1D       - Unweighted
;
;         NOTE: the following special cases apply:
;                * if WEIGHTS is zero, then the corresponding data point
;                  is ignored
;                * if WEIGHTS is NaN or Infinite, and the NAN keyword is
;                  set, then the corresponding data point is ignored
;                * if WEIGHTS is negative, then the absolute value of
;                  WEIGHTS is used.
;
;   XTOL - a nonnegative input variable. Termination occurs when the
;          relative error between two consecutive iterates is at most
;          XTOL (and STATUS is accordingly set to 2 or 3).  Therefore,
;          XTOL measures the relative error desired in the approximate
;          solution.  Default: 1D-10
;
;   YFIT - the best-fit model function, as returned by MYFUNCT.
;
;
; EXAMPLE:
;
;   ; First, generate some synthetic data
;   npts = 200
;   x  = dindgen(npts) * 0.1 - 10.                  ; Independent variable
;   yi = gauss1(x, [2.2D, 1.4, 3000.])              ; "Ideal" Y variable
;   y  = yi + randomn(seed, npts) * sqrt(1000. + yi); Measured, w/ noise
;   sy = sqrt(1000.D + y)                           ; Poisson errors
;
;   ; Now fit a Gaussian to see how well we can recover
;   p0 = [1.D, 1., 1000.]                   ; Initial guess (cent, width, area)
;   p = mpfitfun('GAUSS1', x, y, sy, p0)    ; Fit a function
;   print, p
;
;   Generates a synthetic data set with a Gaussian peak, and Poisson
;   statistical uncertainty.  Then the same function is fitted to the
;   data (with different starting parameters) to see how close we can
;   get.
;
;
; COMMON BLOCKS:
;
;   COMMON MPFIT_ERROR, ERROR_CODE
;
;     User routines may stop the fitting process at any time by
;     setting an error condition.  This condition may be set in either
;     the user's model computation routine (MYFUNCT), or in the
;     iteration procedure (ITERPROC).
;
;     To stop the fitting, the above common block must be declared,
;     and ERROR_CODE must be set to a negative number.  After the user
;     procedure or function returns, MPFIT checks the value of this
;     common block variable and exits immediately if the error
;     condition has been set.  By default the value of ERROR_CODE is
;     zero, indicating a successful function/procedure call.
;
; REFERENCES:
;
;   MINPACK-1, Jorge More', available from netlib (www.netlib.org).
;   "Optimization Software Guide," Jorge More' and Stephen Wright,
;     SIAM, *Frontiers in Applied Mathematics*, Number 14.
;
; MODIFICATION HISTORY:
;   Written, Apr-Jul 1998, CM
;   Added PERROR keyword, 04 Aug 1998, CM
;   Added COVAR keyword, 20 Aug 1998, CM
;   Added ITER output keyword, 05 Oct 1998
;      D.L Windt, Bell Labs, windt@bell-labs.com;
;   Added ability to return model function in YFIT, 09 Nov 1998
;   Analytical derivatives allowed via AUTODERIVATIVE keyword, 09 Nov 1998
;   Parameter values can be tied to others, 09 Nov 1998
;   Cosmetic documentation updates, 16 Apr 1999, CM
;   More cosmetic documentation updates, 14 May 1999, CM
;   Made sure to update STATUS, 25 Sep 1999, CM
;   Added WEIGHTS keyword, 25 Sep 1999, CM
;   Changed from handles to common blocks, 25 Sep 1999, CM
;     - commons seem much cleaner and more logical in this case.
;   Alphabetized documented keywords, 02 Oct 1999, CM
;   Added QUERY keyword and query checking of MPFIT, 29 Oct 1999, CM
;   Corrected EXAMPLE (offset of 1000), 30 Oct 1999, CM
;   Check to be sure that X and Y are present, 02 Nov 1999, CM
;   Documented PERROR for unweighted fits, 03 Nov 1999, CM
;   Changed to ERROR_CODE for error condition, 28 Jan 2000, CM
;   Corrected errors in EXAMPLE, 26 Mar 2000, CM
;   Copying permission terms have been liberalized, 26 Mar 2000, CM
;   Propagated improvements from MPFIT, 17 Dec 2000, CM
;   Added CASH statistic, 10 Jan 2001
;   Added NFREE and NPEGGED keywords, 11 Sep 2002, CM
;   Documented RELSTEP field of PARINFO (!!), CM, 25 Oct 2002
;   Add DOF keyword to return degrees of freedom, CM, 23 June 2003
;   Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
;   Move STRICTARR compile option inside each function/procedure, 9
;     Oct 2006
;   Add NAN keyword, to ignore non-finite data values, 28 Oct 2006, CM
;   Clarify documentation on user-function, derivatives, and PARINFO,
;     27 May 2007
;   Fix bug in handling of explicit derivatives with errors/weights
;     (the weights were not being applied), CM, 03 Sep 2007
;   Add COMPATIBILITY section, CM, 13 Dec 2007
;   Add documentation about NAN behavior, CM, 30 Mar 2009
;
;  $Id: mpfitfun.pro,v 1.15 2009/03/30 16:27:43 craigm Exp $
;-
; Copyright (C) 1997-2002, 2003, 2006, 2007, 2009, Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-

FORWARD_FUNCTION mpfitfun_eval, mpfitfun, mpfit

; This is the call-back function for MPFIT.  It evaluates the
; function, subtracts the data, and returns the residuals.
function mpfitfun_eval, p, dp, _EXTRA=extra

  COMPILE_OPT strictarr
  common mpfitfun_common, fcn, x, y, err, wts, f, fcnargs

  ;; Save the original DP matrix for later use
  if n_params() GT 1 then if n_elements(dp) GT 0 then dp0 = dp

  ;; The function is evaluated here.  There are four choices,
  ;; depending on whether (a) FUNCTARGS was passed to MPFITFUN, which
  ;; is passed to this function as "hf"; or (b) the derivative
  ;; parameter "dp" is passed, meaning that derivatives should be
  ;; calculated analytically by the function itself.
  if n_elements(fcnargs) GT 0 then begin
      if n_params() GT 1 then f = call_function(fcn, x, p, dp, _EXTRA=fcnargs)$
      else                    f = call_function(fcn, x, p, _EXTRA=fcnargs)
  endif else begin
      if n_params() GT 1 then f = call_function(fcn, x, p, dp) $
      else                    f = call_function(fcn, x, p)
  endelse

  np = n_elements(p)
  nf = n_elements(f)

  ;; Compute the deviates, applying either errors or weights
  if n_elements(wts) GT 0 then begin
      result = (y-f)*wts
      if n_elements(dp0) GT 0 AND n_elements(dp) EQ np*nf then begin
          for j = 0L, np-1 do dp[j*nf] = dp[j*nf:j*nf+nf-1] * wts
      endif
  endif else if n_elements(err) GT 0 then begin
      result = (y-f)/err
      if n_elements(dp0) GT 0 AND n_elements(dp) EQ np*nf then begin
          for j = 0L, np-1 do dp[j*nf] = dp[j*nf:j*nf+nf-1] / err
      endif
  endif else begin
      result = (y-f)
  endelse

  ;; Make sure the returned result is one-dimensional.
  result = reform(result, n_elements(result), /overwrite)
  return, result

end

;; Implement residual and gradient scaling according to the
;; prescription of Cash (ApJ, 228, 939)
pro mpfitfun_cash, resid, dresid
  COMPILE_OPT strictarr
  common mpfitfun_common, fcn, x, y, err, wts, f, fcnargs

  sz = size(dresid)
  m = sz[1]
  n = sz[2]

  ;; Do rudimentary dimensions checks, so we don't do something stupid
  if n_elements(y) NE m OR n_elements(f) NE m OR n_elements(resid) NE m then begin
      DIM_ERROR:
      message, 'ERROR: dimensions of Y, F, RESID or DRESID are not consistent'
  endif

  ;; Scale gradient by sqrt(y)/f
  gfact = temporary(dresid) * rebin(reform(sqrt(y)/f,m,1),m,n)
  dresid = reform(dresid, m, n, /overwrite)

  ;; Scale residuals by 1/sqrt(y)
  resid = temporary(resid)/sqrt(y)

  return
end

function mpfitfun, fcn, x, y, err, p, WEIGHTS=wts, FUNCTARGS=fa, $
                   BESTNORM=bestnorm, nfev=nfev, STATUS=status, $
                   parinfo=parinfo, query=query, CASH=cash, $
                   covar=covar, perror=perror, yfit=yfit, $
                   niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                   quiet=quiet, ERRMSG=errmsg, NAN=NAN, _EXTRA=extra

  COMPILE_OPT strictarr
  status = 0L
  errmsg = ''

  ;; Detect MPFIT and crash if it was not found
  catch, catcherror
  if catcherror NE 0 then begin
      MPFIT_NOTFOUND:
      catch, /cancel
      message, 'ERROR: the required function MPFIT must be in your IDL path', /info
      return, !values.d_nan
  endif
  if mpfit(/query) NE 1 then goto, MPFIT_NOTFOUND
  catch, /cancel
  if keyword_set(query) then return, 1

  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = MPFITFUN('MYFUNCT', X, Y, ERR, "+ $
        "START_PARAMS, ... )", /info
      return, !values.d_nan
  endif
  if n_elements(x) EQ 0 OR n_elements(y) EQ 0 then begin
      message, 'ERROR: X and Y must be defined', /info
      return, !values.d_nan
  endif

  if n_elements(err) GT 0 OR n_elements(wts) GT 0 AND keyword_set(cash) then begin
      message, 'ERROR: WEIGHTS or ERROR cannot be specified with CASH', /info
      return, !values.d_nan
  endif
  if keyword_set(cash) then begin
      scalfcn = 'mpfitfun_cash'
  endif

  ;; Use common block to pass data back and forth
  common mpfitfun_common, fc, xc, yc, ec, wc, mc, ac
  fc = fcn & xc = x & yc = y & mc = 0L
  ;; These optional parameters must be undefined first
  ac = 0 & dummy = size(temporary(ac))
  ec = 0 & dummy = size(temporary(ec))
  wc = 0 & dummy = size(temporary(wc))

  ;; FUNCTARGS
  if n_elements(fa) GT 0 then ac = fa

  ;; WEIGHTS or ERROR
  if n_elements(wts) GT 0 then begin
      wc = sqrt(abs(wts))
  endif else if n_elements(err) GT 0 then begin
      wh = where(err EQ 0, ct)
      if ct GT 0 then begin
          errmsg = 'ERROR: ERROR value must not be zero.  Use WEIGHTS instead.'
          message, errmsg, /info
          return, !values.d_nan
      endif
      ;; Appropriate weight for gaussian errors
      wc = 1/abs(err)
  endif

  ;; Check for weights/errors which do not match the dimension
  ;; of the data points
  if n_elements(wc) GT 0 AND $
    n_elements(wc) NE 1 AND $
    n_elements(wc) NE n_elements(yc) then begin
      errmsg = 'ERROR: ERROR/WEIGHTS must either be a scalar or match the number of Y values'
      message, errmsg, /info
      return, !values.d_nan
  endif

  ;; If the weights/errors are a scalar value, and not finite, then
  ;; the fit will surely fail
  if n_elements(wc) EQ 1 then begin
      if finite(wc[0]) EQ 0 then begin
          errmsg = 'ERROR: the supplied scalar WEIGHT/ERROR value was not finite'
          message, errmsg, /info
          return, !values.d_nan
      endif
  endif

  ;; Handle the cases of non-finite data points or weights
  if keyword_set(nan) then begin
      ;; Non-finite data points
      wh = where(finite(yc) EQ 0, ct)
      if ct GT 0 then begin
          yc[wh] = 0
          ;; Careful: handle case when weights were a scalar...
          ;;   ... promote to a vector
          if n_elements(wc) EQ 1 then wc = replicate(wc[0], n_elements(yc))
          wc[wh] = 0
      endif

      ;; Non-finite weights
      wh = where(finite(wc) EQ 0, ct)
      if ct GT 0 then wc[wh] = 0
  endif

  result = mpfit('mpfitfun_eval', p, SCALE_FCN=scalfcn, $
                 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                 covar=covar, perror=perror, $
                 niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)

  ;; Retrieve the fit value
  yfit = temporary(mc)
  ;; Some cleanup
  xc = 0 & yc = 0 & wc = 0 & ec = 0 & mc = 0 & ac = 0

  ;; Print error message if there is one.
  if NOT keyword_set(quiet) AND errmsg NE '' then $
    message, errmsg, /info

  return, result
end


pro sxdelpar, h, parname
;+
; NAME:
;   SXDELPAR
; PURPOSE:
;   Procedure to delete a keyword parameter(s) from a FITS header
;
; CALLING SEQUENCE:
;   sxdelpar, h, parname
;
; INPUTS:
;   h - FITS or STSDAS header, string array
;   parname - string or string array of keyword name(s) to delete
;
; OUTPUTS:
;   h - updated FITS header, If all lines are deleted from
;     the header, then h is returned with a value of 0
;
; EXAMPLE:
;   Delete the astrometry keywords CDn_n from a FITS header, h
;
;   IDL> sxdelpar, h, ['CD1_1','CD1_2','CD2_1','CD2_2']
;
; NOTES:
;   (1)  No message is returned if the keyword to be deleted is not found
;   (2)  All appearances of a keyword in the header will be deleted
; HISTORY:
;   version 1  D. Lindler Feb. 1987
;   Converted to new IDL  April 1990 by D. Lindler
;   Test for case where all keywords are deleted    W. Landsman Aug 1995
;   Converted to IDL V5.0   W. Landsman   September 1997
;       Allow for headers with more than 32767 lines W. Landsman Jan. 2003
;------------------------------------------------------------------
 On_error,2

 if N_Params() LT 2 then begin
      print,'Syntax - SXDELPAR, h, parname'
      return
 endif

; convert parameters to string array of upper case names of length 8 char

 s = size(parname) & ndim = s[0] & type = s[ndim+1]

 if type NE 7 then message,'Keyword name(s) must be a string or string array'

 num = N_elements( parname )
 par = strtrim( strupcase(parname),2 )

 s = size(h) & ndim = s[0] & type = s[ndim+1]

 if (ndim NE 1) or (type NE 7) then $
    message,'FITS header (1st parameter) must be a string array'

 nlines = s[1]   ;number of lines in header array
 pos = 0L     ;position in compressed header with keywords removed

; loop on header lines

 keyword = strtrim( strmid(h,0,8), 2 )
 for i = 0L, nlines-1 do begin
    for j = 0,num-1 do $
              if keyword[i] EQ par[j] then goto, DELETE ;delete it?
    h[pos] = h[i]   ;keep it
    pos = pos+1       ;increment number of lines kept
    if keyword[i] eq 'END' then goto, DONE     ;end of header
DELETE:
 endfor

DONE:
 if pos GT 0 then h = h[0:pos-1] else h = 0       ;truncate

 return
 end



FUNCTION CONTINUUM, LAMBDA, SPLINE_TENSION, GRAPH=GRAPH

; purpose:
; calculate solar spectrum continuum for given wavelenght
;
; usage:
; result = continuum (wavelenght [, spline_tension, /GRAPH])

; input:
; wavelenght in angstoems

; output:
; intensity of solar continuum in  egr * sr^(-1) * cm^(-2) * A^(-1) * s^(-1)

; references:
; Allen

; author
; Maciek Zapior

;table from Allen p.172 column 5:
I = [   0.04, 0.20,   0.30, 0.5,    1.19,  2.15,    2.83,  3.01,    3.20,  $
       3.62, 4.1,    4.4,   4.58, 4.60,   4.59, 4.55,   4.54, 4.48,   $
       4.40, 4.31,   4.08, 3.68,   3.27, 2.88,   2.53, 2.24,   $
       1.97, 1.55,   1.23, 0.99,   0.81, 0.564,  0.403,   0.268,    $
       0.183,    0.081, 0.041,  0.0135,  0.0057   ]

;table from Allen p.172 column 1:
MI =[   0.20, 0.22,   0.24, 0.26,   0.28, 0.30,   0.32, 0.34,   0.36, $
       0.37, 0.38,   0.39, 0.40,   0.41, 0.42,   0.43, 0.44,   0.45, $
       0.46, 0.48,   0.50, 0.55,   0.60, 0.65,   0.70, 0.75,   $
       0.80, 0.9,    1.0,   1.1,  1.2, 1.4,    1.6,   1.8,  $
       2.0,  2.5, 3.0,    4.0,   5.0     ]

LAMBDA = ROUND(LAMBDA)

IF LAMBDA LT 2000 OR LAMBDA GT 50000 THEN BEGIN
    PRINT, 'INPUT VALUE MUST BE BETWEEN 2000 AND 50000'
    STOP
ENDIF

if N_PARAMS() EQ 1 THEN BEGIN
    SPLINE_TENSION = 10.0
    ENDIF
XX = FINDGEN(50001)/10000.

I_SPL = SPLINE(MI, I, XX, SPLINE_TENSION)

IF KEYWORD_SET(GRAPH) THEN BEGIN
    LOADCT,38
    WINDOW, 0, XSIZE=1280, YSIZE=640
    PLOT, MI, I, PSYM=2, XRANGE = [0.2,0.9], XSTYLE=1
    OPLOT, XX, I_SPL, COLOR=150
    PLOTS, [LAMBDA/10000.,LAMBDA/10000.],[I_SPL(LAMBDA)-0.5, I_SPL(LAMBDA)+0.5], COLOR=200
    WSHOW, 0

ENDIF

RETURN, I_SPL(LAMBDA) * 10.^6

END


FUNCTION DARKENING, THETA, LAMBDA, SPLINE_TENSION, GRAPH=GRAPH

;   purpose:
; calculate limb darkening of the Sun using Allen's formula
;
;    usage:
; RESULT = DARKENING ( theta(radians), wavelenght(angstroems) [, spline tension, /GRAPH] )
;

; output:
; factor of limb darkening for given wavelenght in respect to disc center

; References:
; Allen
;
; author:
; Maciek Zapior

if N_PARAMS() EQ 2 THEN SPLINE_TENSION = 1.0

; table taken from Allen:

MI= [   0.20, 0.22,   0.245,   0.265,    0.28,       $
       0.30, 0.32,   0.35, 0.37,   0.38,        $
       0.40, 0.45,   0.50, 0.55,   0.60, 0.80,   $
       1.0,  1.5, 2.0,    3.0,   5.0,  10.0 ]

U = [   0.12, -1.3,   -0.1,    -0.1,  0.38,       $
        0.74, 0.88,   0.98, 1.03,   0.92,          $
        0.91, 0.99,   0.97,    0.93,  0.88,    0.73,  $
        0.64, 0.57,   0.48, 0.35,   0.22, 0.15    ]

V = [   0.33, 1.6,    0.85,  0.90,    0.57,         $
       0.20, 0.03,   -0.10,    -0.16, -0.05,       $
       -0.05,    -0.17, -0.22,  -0.23,   -0.23,    -0.22,     $
       -0.20,    -0.21, -0.18,  -0.12,   -0.07,    -0.07 ]

;creating abscissa vector
XX = INDGEN(10000) ; ABSCISSA
LAMBDA_TABLE = MI * 10.0^4 ; RECALCULATING WAVELENGHT TO ANGSTROMS

;spline interpolation for U and V coefficients
SPL_V = SPLINE(LAMBDA_TABLE, V, XX, SPLINE_TENSION)
SPL_U = SPLINE(LAMBDA_TABLE, U, XX, SPLINE_TENSION)

;calculation of limb darkening factor using Allen's formula for wavelenght eq LAMBDA and angle eq THETA
DARKENING_FACTOR = 1.0 - SPL_U(LAMBDA) - SPL_V(LAMBDA) + SPL_U(LAMBDA)*COS(THETA) + SPL_V(LAMBDA)*(COS(THETA))^2.0

IF KEYWORD_SET(GRAPH) THEN BEGIN

    WINDOW, 0, XSIZE = 1280, YSIZE=400, XPOS=0, YPOS=0, TITLE='V2 COEFFICIENT'
    SPL_V = SPLINE(LAMBDA_TABLE, V, XX, SPLINE_TENSION)
    PLOT, LAMBDA_TABLE, V, XRANGE=[3800, 9000], YSTYLE=1, YRANGE=[-0.3, 0.3],PSYM=2
    OPLOT, XX, SPL_V, COLOR = 150
    WSHOW, 0


    WINDOW, 1, XSIZE = 1280, YSIZE=400, XPOS=0, YPOS=400, TITLE='U2 COEFFICIENT'
    SPL_U = SPLINE(LAMBDA_TABLE, U, XX, SPLINE_TENSION)
    PLOT, LAMBDA_TABLE, U, XRANGE=[3800, 9000], YSTYLE=1, YRANGE=[0.5, 1.2], PSYM=2
    OPLOT, XX, SPL_U, COLOR = 150
    WSHOW, 1

ENDIF

RETURN, DARKENING_FACTOR

END




;+
; NAME:
;   MPFITELLIPSE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE:
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Approximate fit to points forming an ellipse
;
; MAJOR TOPICS:
;   Curve and Surface Fitting
;
; CALLING SEQUENCE:
;   parms = MPFITELLIPSE(X, Y, start_parms, [/TILT, WEIGHTS=wts, ...])
;
; DESCRIPTION:
;
;   MPFITELLIPSE fits a closed elliptical or circular curve to a two
;   dimensional set of data points.  The user specifies the X and Y
;   positions of the points, and an optional set of weights.  The
;   ellipse may also be tilted at an arbitrary angle.
;
;   IMPORTANT NOTE: this fitting program performs simple ellipse
;   fitting.  It will not work well for ellipse data with high
;   eccentricity.  More robust answers can usually be obtained with
;   "orthogonal distance regression."  (See FORTRAN package ODRPACK on
;   netlib.org for more information).
;
;   The best fitting ellipse parameters are returned from by
;   MPFITELLIPSE as a vector, whose values are:
;
;      P[0]   Ellipse semi axis 1
;      P[1]   Ellipse semi axis 2   ( = P[0] if CIRCLE keyword set)
;      P[2]   Ellipse center - x value
;      P[3]   Ellipse center - y value
;      P[4]   Ellipse rotation angle (radians) if TILT keyword set
;
;   As of the 17 Sep 2007 version of this function, P[0] is meant to
;   be the semi-major axis, and P[1] the semi-minor axis.  The
;   returned semi-axis lengths should always be positive.
;
;   The user may specify an initial set of trial parameters, but by
;   default MPFITELLIPSE will estimate the parameters automatically.
;
;   Users should be aware that in the presence of large amounts of
;   noise, namely when the measurement error becomes significant
;   compared to the ellipse axis length, then the estimated parameters
;   become unreliable.  Generally speaking the computed axes will
;   overestimate the true axes.  For example when (SIGMA_R/R) becomes
;   0.5, the radius of the ellipse is overestimated by about 40%.
;
;   This unreliability is also pronounced if the ellipse has high
;   eccentricity, as noted above.
;
;   Users can weight their data as they see appropriate.  However, the
;   following prescription for the weighting may serve as a good
;   starting point, and appeared to produce results comparable to the
;   typical chi-squared value.
;
;     WEIGHTS = 0.75/(SIGMA_X^2 + SIGMA_Y^2)
;
;   where SIGMA_X and SIGMA_Y are the measurement error vectors in the
;   X and Y directions respectively.  However, this has not been
;   robustly tested, and it should be pointed out that this weighting
;   may only be appropriate for a set of points whose measurement
;   errors are comparable.  If a more robust estimation of the
;   parameter values is needed, the so-called orthogonal distance
;   regression package should be used (ODRPACK, available in FORTRAN
;   at www.netlib.org).
;
; INPUTS:
;
;   X - measured X positions of the points in the ellipse.
;   Y - measured Y positions of the points in the ellipse.
;
;   START_PARAMS - an array of starting values for the ellipse
;                  parameters, as described above.  This parameter is
;                  optional; if not specified by the user, then the
;                  ellipse parameters are estimated automatically from
;                  the properties of the data.
;
; RETURNS:
;
;   Returns the best fitting model ellipse parameters.
;
; KEYWORDS:
;
;   ** NOTE ** Additional keywords such as PARINFO, BESTNORM, and
;              STATUS are accepted by MPFITELLIPSE but not documented
;              here.  Please see the documentation for MPFIT for the
;              description of these advanced options.
;
;   PERROR - upon return, the 1-sigma uncertainties of the returned
;            ellipse parameter values.  These values are only
;            meaningful if the WEIGHTS keyword is specified properly.
;
;            If the fit is unweighted (i.e. no errors were given, or
;            the weights were uniformly set to unity), then PERROR
;            will probably not represent the true parameter
;            uncertainties.
;
;   QUIET - if set then diagnostic fitting messages are suppressed.
;           Default: QUIET=1 (i.e., no diagnostics]
;
;   CIRCULAR - if set, then the curve is assumed to be a circle
;              instead of ellipse.  When set, the parameters P[0] and
;              P[1] will be identical and the TILT keyword will have
;              no effect.
;
;   TILT - if set, then the major and minor axes of the ellipse
;          are allowed to rotate with respect to the data axes.
;          Parameter P[4] will be set to the clockwise rotation angle
;          of the P[0] axis in radians, as measured from the +X axis.
;          P[4] should be in the range 0 to !dpi.
;
;   WEIGHTS - Array of weights to be used in calculating the
;             chi-squared value.  The chi-squared value is computed
;             as follows:
;
;                CHISQ = TOTAL( (Z-MYFUNCT(X,Y,P))^2 * ABS(WEIGHTS)^2 )
;
;             Users may wish to follow the guidelines for WEIGHTS
;             described above.
;
;
; EXAMPLE:
;
; ; Construct a set of points on an ellipse, with some noise
;   ph0 = 2*!pi*randomu(seed,50)
;   x =  50. + 32.*cos(ph0) + 4.0*randomn(seed, 50)
;   y = -75. + 65.*sin(ph0) + 0.1*randomn(seed, 50)
;
; ; Compute weights function
;   weights = 0.75/(4.0^2 + 0.1^2)
;
; ; Fit ellipse and plot result
;   p = mpfitellipse(x, y)
;   plot, x, y, psym=1
;   phi = dindgen(101)*2D*!dpi/100
;   oplot, p[2]+p[0]*cos(phi), p[3]+p[1]*sin(phi)
;
; REFERENCES:
;
;   MINPACK-1, Jorge More', available from netlib (www.netlib.org).
;   "Optimization Software Guide," Jorge More' and Stephen Wright,
;     SIAM, *Frontiers in Applied Mathematics*, Number 14.
;
; MODIFICATION HISTORY:
;
;   Ported from MPFIT2DPEAK, 17 Dec 2000, CM
;   More documentation, 11 Jan 2001, CM
;   Example corrected, 18 Nov 2001, CM
;   Change CIRCLE keyword to the correct CIRCULAR keyword, 13 Sep
;      2002, CM
;   Add error messages for SYMMETRIC and CIRCLE, 08 Nov 2002, CM
;   Found small error in computation of _EVAL (when CIRCULAR) was set;
;      sanity check when CIRCULAR is set, 21 Jan 2003, CM
;   Convert to IDL 5 array syntax (!), 16 Jul 2006, CM
;   Move STRICTARR compile option inside each function/procedure, 9
;     Oct 2006
;   Add disclaimer about the suitability of this program for fitting
;     ellipses, 17 Sep 2007, CM
;   Clarify documentation of TILT angle; make sure output contains
;    semi-major axis first, followed by semi-minor; make sure that
;    semi-axes are always positive (and can handle negative inputs)
;      17 Sep 2007, CM
;   Output tilt angle is now in range 0 to !DPI, 20 Sep 2007, CM
;   Some documentation clarifications, including to remove reference
;     to the "ERR" keyword, which does not exist, 17 Jan 2008, CM
;
;  $Id: mpfitellipse.pro,v 1.12 2008/01/17 21:18:36 craigm Exp $
;-
; Copyright (C) 1997-2000,2002,2003,2007,2008 Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy, modify, and distribute modified or
; unmodified copies is granted, provided this copyright and disclaimer
; are included unchanged.
;-



PRO C_D, R

WINDOW, 31, XSIZE=2000, YSIZE=1500,xpos=-50, ypos=-50,TITLE=''

FOR I=0,30 DO BEGIN
       WINDOW,I,XSIZE=RANDOMU(SEED)*1200+100,YSIZE=RANDOMU(SEED1)*800+100,XPOS=RANDOMU(SEED2)*500,YPOS=RANDOMU(SEED3)*500,TITLE=''
       WAIT,0.2
ENDFOR

MENU,['YES','NO'],$
    SO,TITLE='IDL FOUND FATAL ERROR, DO YOU WANT EXIT THE PROGRAM?',BUTSIZE=600,xoffset=450,yoffset=350

    WINDOW, 31, XSIZE=2000, YSIZE=1500,xpos=-50, ypos=-50,TITLE=''
    RRR = bytarr(2000,1500)
    FOR I = 0, 50 DO BEGIN
    RRR(*,*) = RRR(*,*) + randomu(seed)*10.
    TV, RRR
    loadct,i mod 40
ENDFOR

    MENU,['OK','CANCEL'],$
    SO,TITLE='EXIT IDL?',BUTSIZE=600,xoffset=450,yoffset=350

REPEAT BEGIN
    FOR I=0,30 DO BEGIN
       WINDOW,I,XSIZE=RANDOMU(SEED)*1200+100,YSIZE=RANDOMU(SEED1)*800+100,XPOS=RANDOMU(SEED2)*500,YPOS=RANDOMU(SEED3)*500,TITLE=''
       RR = DIST(1200,1200)
       RR = RR * RANDOMU(SEED5)*255
       LOADCT,round(randomu(seed7)*40)
       TVSCL,RR
       WAIT,0.4*RANDOMU(SEED6)
    ENDFOR

    MENU,['YES','NO','I DONT KNOW','HELP ME PLEASE','I AM REALLY SORRY I PRESSED THIS BUTTON'],$
    SORRY,TITLE='      I SAID DONT PRESS THIS BUTTON !!!',BUTSIZE=500,xoffset=450,yoffset=350

ENDREP UNTIL SORRY EQ 4 OR SORRY EQ 5

WINDOW,1,XSIZE=500,YSIZE=0,XPOS=500,YPOS=100,TITLE='                BEST REGARDS FROM MACIEJ ZAPIOR    ;) '

WAIT,2

FOR I=0,31 DO WDELETE,I

LOADCT,3

END


FORWARD_FUNCTION mpfitellipse_u, mpfitellipse_eval, mpfitellipse, mpfit

; Compute the "u" value = (x/a)^2 + (y/b)^2 with optional rotation
function mpfitellipse_u, x, y, p, tilt=tilt, circle=circle
  COMPILE_OPT strictarr
  widx  = abs(p[0]) > 1e-20 & widy  = abs(p[1]) > 1e-20
  if keyword_set(circle) then widy  = widx
  xp    = x-p[2]            & yp    = y-p[3]
  theta = p[4]

  if keyword_set(tilt) AND theta NE 0 then begin
      c  = cos(theta) & s  = sin(theta)
      return, ( (xp * (c/widx) - yp * (s/widx))^2 + $
                (xp * (s/widy) + yp * (c/widy))^2 )
  endif else begin
      return, (xp/widx)^2 + (yp/widy)^2
  endelse

end

; This is the call-back function for MPFIT.  It evaluates the
; function, subtracts the data, and returns the residuals.
function mpfitellipse_eval, p, tilt=tilt, circle=circle, _EXTRA=extra

  COMPILE_OPT strictarr
  common mpfitellipse_common, xy, wc

  tilt = keyword_set(tilt)
  circle = keyword_set(circle)
  u2 = mpfitellipse_u(xy[*,0], xy[*,1], p, tilt=tilt, circle=circle) - 1.

  if n_elements(wc) GT 0 then begin
      if circle then u2 = sqrt(abs(p[0]*p[0]*wc))*u2 $
      else           u2 = sqrt(abs(p[0]*p[1]*wc))*u2
  endif

  return, u2
end

function mpfitellipse, x, y, p0, WEIGHTS=wts, $
                       BESTNORM=bestnorm, nfev=nfev, STATUS=status, $
                       tilt=tilt, circular=circle, $
                       circle=badcircle1, symmetric=badcircle2, $
                       parinfo=parinfo, query=query, $
                       covar=covar, perror=perror, niter=iter, $
                       quiet=quiet, ERRMSG=errmsg, _EXTRA=extra

  COMPILE_OPT strictarr
  status = 0L
  errmsg = ''

  ;; Detect MPFIT and crash if it was not found
  catch, catcherror
  if catcherror NE 0 then begin
      MPFIT_NOTFOUND:
      catch, /cancel
      message, 'ERROR: the required function MPFIT must be in your IDL path', /info
      return, !values.d_nan
  endif
  if mpfit(/query) NE 1 then goto, MPFIT_NOTFOUND
  catch, /cancel
  if keyword_set(query) then return, 1

  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = MPFITELLIPSE(X, Y, START_PARAMS, ... )", $
        /info
      return, !values.d_nan
  endif
  nx = n_elements(x) & ny = n_elements(y)
  if (nx EQ 0) OR (ny EQ 0) OR (nx NE ny) then begin
      message, 'ERROR: X and Y must have the same number of elements', /info
      return, !values.d_nan
  endif

  if keyword_set(badcircle1) OR keyword_set(badcircle2) then $
    message, 'ERROR: do not use the CIRCLE or SYMMETRIC keywords.  ' +$
    'Use CIRCULAR instead.'

  p = make_array(5, value=x[0]*0)

  if n_elements(p0) GT 0 then begin
      p[0] = p0
      if keyword_set(circle) then p[1] = p[0]
  endif else begin
      mx = moment(x)
      my = moment(y)
      p[0] = [sqrt(mx[1]), sqrt(my[1]), mx[0], my[0], 0]
      if keyword_set(circle) then $
        p[0:1] = sqrt(mx[1]+my[1])
  endelse

  common mpfitellipse_common, xy, wc
  if n_elements(wts) GT 0 then begin
      wc = abs(wts)
  endif else begin
      wc = 0 & dummy = temporary(wc)
  endelse

  xy = [[x],[y]]

  result = mpfit('mpfitellipse_eval', p, $
                 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                 covar=covar, perror=perror, niter=iter, $
                 functargs={circle:keyword_set(circle), tilt:keyword_set(tilt)},$
                 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)

  ;; Print error message if there is one.
  if NOT keyword_set(quiet) AND errmsg NE '' then $
    message, errmsg, /info

  ;; Sanity check on resulting parameters
  if keyword_set(circle) then begin
      result[1] = result[0]
      perror[1] = perror[0]
  endif
  if NOT keyword_set(tilt) then begin
      result[4] = 0
      perror[4] = 0
  endif

  ;; Make sure the axis lengths are positive, and the semi-major axis
  ;; is listed first
  result[0:1] = abs(result[0:1])
  if abs(result[0]) LT abs(result[1]) then begin
      tmp = result[0] & result[0] = result[1] & result[1] = tmp
      tmp = perror[0] & perror[0] = perror[1] & perror[1] = tmp
      if keyword_set(tilt) then result[4] = result[4] - !dpi/2d
  endif

  if keyword_set(tilt) then begin
      ;; Put tilt in the range 0 to +Pi
      result[4] = result[4] - !dpi * floor(result[4]/!dpi)
  endif

  return, result
end


;----------------    FITTENG CURVE TO DATA
FORWARD_FUNCTION mpevalexpr, mpfitexpr_eval, mpfitexpr, mpfit

; Utility function which simply returns the value of the expression,
; evaluated at each point in x, using the parameters p.
function mpevalexpr, _expr, x, p, functargs=private

  COMPILE_OPT strictarr
  _cmd = '_f = '+_expr
  _err = execute(_cmd)
  return, _f
end

; This is the call-back function for MPFIT.  It evaluates the
; expression, subtracts the data, and returns the residuals.
function mpfitexpr_eval, p, _EXTRA=private

  COMPILE_OPT strictarr
  common mpfitexpr_common, _expr, x, y, err, _wts, _f

  ;; Compute the model value by executing the expression
  _f = 0.D
  _cmd = '_f = '+_expr
  _xxx = execute(_cmd)
  if _xxx EQ 0 then message, 'ERROR: command execution failed.'

  ;; Compute the deviates, applying either errors or weights
  if n_elements(err) GT 0 then begin
      result = (y-_f)/err
  endif else if n_elements(_wts) GT 0 then begin
      result = (y-_f)*_wts
  endif else begin
      result = (y-_f)
  endelse

  ;; The returned result should be one-dimensional
  result = reform(result, n_elements(result), /overwrite)
  return, result
end

;; This is the main entry point for this module
function mpfitexpr, expr, x, y, err, p, WEIGHTS=wts, $
                    BESTNORM=bestnorm, STATUS=status, nfev=nfev, $
                    parinfo=parinfo, query=query, functargs=fcnargs, $
                    covar=covar, perror=perror, yfit=yfit, $
                    niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                    quiet=quiet, _EXTRA=extra, errmsg=errmsg

  COMPILE_OPT strictarr
  status = 0L
  errmsg = ''

  ;; Detect MPFIT and crash if it was not found
  catch, catcherror
  if catcherror NE 0 then begin
      MPFIT_NOTFOUND:
      catch, /cancel
      message, 'ERROR: the required function MPFIT must be in your IDL path', /info
      return, !values.d_nan
  endif
  if mpfit(/query) NE 1 then goto, MPFIT_NOTFOUND
  catch, /cancel
  if keyword_set(query) then return, 1

  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = MPFITEXPR('EXPR', X, Y, ERR, "+ $
        "START_PARAMS, ... )", /info
      return, !values.d_nan
  endif
  if n_elements(x) EQ 0 OR n_elements(y) EQ 0 then begin
      message, 'ERROR: X and Y must be defined', /info
      return, !values.d_nan
  endif

  ;; If no parameters are given, then parse the input expression,
  ;; and determine the number of parameters automatically.
  if (n_elements(parinfo) GT 0) AND (n_elements(p) EQ 0) then $
    p = parinfo[*].value
  if (n_elements(p) EQ 0) then begin
      pos = 0L
      nparams = 0L
      ee = strupcase(expr)
      ;; These are character constants representing the boundaries of
      ;; variable names.
      ca = (byte('A'))[0]
      cz = (byte('Z'))[0]
      c0 = (byte('0'))[0]
      c9 = (byte('9'))[0]
      c_ = (byte('_'))[0]  ;; Underscore can be in a variable name
      ll = strlen(ee)
      pnames = ['']

      ;; Now step through, looking for variables looking like p[0], etc.
      repeat begin
          i = [strpos(ee, 'P(', pos), strpos(ee, 'P[', pos)]
          wh = where(i GE 0, ct)
          if ct LE 0 then goto, DONE_PARAMS
          i = min(i[wh])

          ;; None found, finished
          if i LT 0 then goto, DONE_PARAMS
          ;; Too close to the end of the string
          if i GT ll-4 then goto, DONE_PARAMS

          ;; Have to be careful here, to be sure that this isn't just
          ;; a variable name ending in "p"
          maybe = 0
          ;; If this is the first character
          if i EQ 0 then maybe = 1 $
          else begin
              ;; Or if the preceding character is a non-variable character
              c = (byte(strmid(ee, i-1, 1)))[0]
              if NOT ( (c GE ca AND c LE cz) OR (c GE c0 AND c LE c9) $
                       OR c EQ c_ ) then maybe = 1
          endelse
          if maybe then begin
              ;; If we found one, then strip out the value inside the
              ;; parentheses.
              rest = strmid(ee, i+2, ll-i-2)
              next = str_sep(rest, ')', /trim)
              next = next[0]
              pnames = [pnames, next]
          endif
          pos = i+1
      endrep until pos GE ll

      DONE_PARAMS:
      if n_elements(pnames) EQ 1 then begin
          message, 'ERROR: no parameters to fit', /info
          return, !values.d_nan
      endif

      ;; Finally, we take the maximum parameter number
      pnames = pnames[1:*]
      nparams = max(long(pnames)) + 1
      if NOT keyword_set(quiet) then $
        message, '  Number of parameters: '+strtrim(nparams,2) $
        + ' (initialized to zero)', /info

      ;; Create a parameter vector, starting at zero
      p = dblarr(nparams)
  endif

  ;; Use common block to pass data back and forth
  common mpfitexpr_common, fc, xc, yc, ec, wc, mc
  fc = expr & xc = x & yc = y & mc = 0L
  ;; These optional parameters must be undefined first
  ec = 0 & dummy = size(temporary(ec))
  wc = 0 & dummy = size(temporary(wc))

  if n_elements(wts) GT 0 then begin
      wc = sqrt(abs(wts))
  endif else if n_elements(err) GT 0 then begin
      wh = where(err EQ 0, ct)
      if ct GT 0 then begin
          message, 'ERROR: ERROR value must not be zero.  Use WEIGHTS.', $
            /info
          return, !values.d_nan
      endif
      ec = err
  endif

  ;; Test out the function, as mpfit would call it, to see if it works
  ;; okay.  There is no sense in calling the fitter if the function
  ;; itself doesn't work.
  catch, catcherror
  if catcherror NE 0 then begin
      CATCH_ERROR:
      catch, /cancel
      message, 'ERROR: execution of "'+expr+'" failed.', /info
      message, '       check syntax and parameter usage', /info
      xc = 0 & yc = 0 & ec = 0 & wc = 0 & ac = 0
      return, !values.d_nan
  endif

  ;; Initialize.  Function that is actually called is mpfitexpr_eval,
  ;; which is a wrapper that sets up the expression evaluation.
  fcn = 'mpfitexpr_eval'

  ;; FCNARGS are passed by MPFIT directly to MPFITEXPR_EVAL.  These
  ;; actually contain the data, the expression, and a slot to return
  ;; the model function.
  fvec = call_function(fcn, p, _EXTRA=fcnargs)
  if n_elements(fvec) EQ 1 then $
    if NOT finite(fvec[0]) then goto, CATCH_ERROR
  ;; No errors caught if reached this stage
  catch, /cancel

  ;; Call MPFIT
  result = mpfit(fcn, p, $
                 parinfo=parinfo, STATUS=status, nfev=nfev, BESTNORM=bestnorm,$
                 covar=covar, perror=perror, functargs=fcnargs, $
                 niter=niter, nfree=nfree, npegged=npegged, dof=dof, $
                 ERRMSG=errmsg, quiet=quiet, _EXTRA=extra)

  ;; Retrieve the fit value
  yfit = temporary(mc)
  ;; Some cleanup
  xc = 0 & yc = 0 & wc = 0 & ec = 0 & mc = 0 & ac = 0

  ;; Print error message if there is one.
  if NOT keyword_set(quiet) AND errmsg NE '' then $
    message, errmsg, /info

  return, result
end


;----------------    END FITTENG CURVE TO DATA



PRO VELOC_MODIFY,LP,PR,V,VMEAN,PL,NO_CAM,LINE_INDEX ;CO TU WCHODZI?

;LP TO TABLICA Z LAMBDAMI (25 PKTOW)
;PR TO TABLICA SYGNALOW, JUZ W 25 PKTACH (INTERPOLOWANA Z 9 KANALOW)
;V TABLICA 4 ELEM DO ROZNYCH POMIAROW PREDKOSCI
;PL EQ 0 WTEDY SIE NIE WYSWIETLA OKNO
;NO_CAM INDEX OKNA Z WYKRESEM WIDMA




CASE LINE_INDEX OF
    1: BEGIN
       L0=6562.808   ;ha
    END
    2: BEGIN
       ;L0=3968.492  ;caII(H) sprawdzic
       L0=3933.682 ; ca II (K)
    END
    3: BEGIN
       L0=5875.618   ;d3 sprawdzic
    END
    4: BEGIN
       L0=4861.342   ;hb
    END
    5: BEGIN
       L0=8542.144   ;ir wstawic
    END

ENDCASE


C=299792.458
IMIN=MIN(PR) ;MAX SYGNAL
IMAX=MAX(PR) ;MIN SYGNAL
DI=(IMAX-IMIN)/12.0  ; DI - ROZNICA MIEDZY MAX I MIN SYGNALEM
DL=LP(1)-LP(0)

LSS=FLTARR(4)
N=N_ELEMENTS (PR)
npom=0
FOR L=0,3 DO BEGIN
  I=IMAX-(L+1)*DI   ; TU ZMIANA Z MIN NA MAX, ZEBY OD GORY LICZYL
  V(L)=0.0
  L1=-1.0
  L2=-1.0
  FOR K=0,N-2 DO BEGIN
      R1=PR(K)-I
      R2=PR(K+1)-I
      RC=PR(K+1)-PR(K)
      IF ((R1 GT 0.0) AND (R2 LT 0.0)) THEN L1=LP(K)-DL*R1/RC
      IF ((R1 LT 0.0) AND (R2 GT 0.0)) THEN L2=LP(K)-DL*R1/RC
   ENDFOR
   IF (L1*L2  GT 0.0) THEN BEGIN
     LS=((L1<L2)+(L1>L2))/2.0
     LSS(L)=LS
     V(L)=C*((LS-L0)/L0)
     NPOM=NPOM+1
   ENDIF
   VMEAN=(V(0)+V(1)+V(2)+V(3))/FLOAT(NPOM)
ENDFOR
IF (PL NE 0) THEN BEGIN
;CANCELLING SHOWING NEW WINDOW
;PLOT WILL BE ON RECENTLY OPENED WINDOW
;window,no_cam+1,ypos=0,xpos=80,xsize=320*4,ysize=500,title=' BISECTIES '
 plot,lp,PR,XSTYLE=1,TITLE='V='+STRCOMPRESS(STRING(VMEAN),/REMOVE_ALL)+' KM/S',CHARSIZE=1.8
FOR L=0,3 DO BEGIN ;BYLO 0,3
  I=IMAX-(L+1)*DI ; TU TEZ ZMIANA ANALOGICZNA Z MIN NA MAX I Z (+) NA (-)
  oplot,[min(lp),max(lp)],[I,I] , color=150        ;POZIOME KRESKI
  oplot,[LSS(L),LSS(L)],[0.9*I,1.1*I], color=150    ; PIONOWE KRESECZKI
 ENDFOR
ENDIF

END


FUNCTION CIRCLE_BIN, xcenter, ycenter, radius
   points = (2 * !PI / 999.0) * FINDGEN(1000)
   x = xcenter + radius * COS(points )
   y = (ycenter + radius * SIN(points ))/2.
   RETURN, TRANSPOSE([[x],[y]])
END



FUNCTION NECKEL, THETA, LAMBDA, SPLINE_TENSION, GRAPH=GRAPH

;   purpose:
; calculate limb darkening of the Sun using Neckel's formula
;
;    usage:
; RESULT = NECKEL ( theta(radians), wavelenght(angstroems) [, spline tension, /GRAPH] )
;

; output:
; factor of limb darkening for given wavelenght in respect to disc center

; References:
; Neckel 1994
;
; author:
; Maciek Zapior

IF (N_PARAMS() EQ 0) OR (N_PARAMS() EQ 1)  THEN BEGIN
    PRINT, 'USAGE:  RESULT = NECKEL ( theta ( radians ),   wavelenght ( angstroems )   [, spline tension, /GRAPH] )'
    STOP
ENDIF

if N_PARAMS() EQ 2 THEN SPLINE_TENSION = 1.0

; table taken from Neckler and Labs SoPh 153, 91 (1994):


NM = [  303.327, 310.843, 320.468, 329.897, 349.947, 365.875, 374.086, 390.915, 401.970, 416.319, $
       427.930, 443.885, 445.125, 457.345, 477.427, 492.905, 519.930, 541.760, 559.950, 579.880, $
       610.975, 640.970, 669.400, 700.875, 748.710, 811.760, 869.600, 948.850, 1046.600,1098.950 ]

A0 = [  0.08011, 0.08160, 0.08833, 0.09188, 0.11012, 0.12828, 0.12579, 0.12995, 0.12323, 0.12814, $
       0.14249, 0.16220, 0.15248, 0.16604, 0.19571, 0.20924, 0.23695, 0.26073, 0.26892, 0.28392, $
       0.30854, 0.33644, 0.34685, 0.37885, 0.40627, 0.42977, 0.45396, 0.47855, 0.49870, 0.51149  ]

A1 = [  0.70695, 0.71609, 0.77285, 0.92459, 1.02168, 1.04969, 0.85402, 0.91836, 1.08648, 1.19947, $
       1.28796, 1.24893, 1.38517, 1.38544, 1.30551, 1.30798, 1.29927, 1.27428, 1.34319, 1.36896, $
       1.36620, 1.30590, 1.37539, 1.25553, 1.22842, 1.25182, 1.25101, 1.19813, 1.21429, 1.19354  ]

A2 = [     0.49910,  0.69685,  0.65382,  0.19604,  -0.10924, 0.17482,  0.54601,  -0.07566, -0.43974, -0.84407,  $
              -1.19564, -0.92165, -1.49615, -1.52275, -1.25845, -1.20411, -1.28034,  -1.30352, -1.58427, -1.75998,  $
          -1.83572, -1.79238, -2.04425, -1.70908, -1.67877, -1.85164, -2.02958,  -1.86296, -2.06976, -2.00174  ]

A3 = [  -0.31080,  -0.87703,  -1.04647,  -0.39546,  -0.00055,  -1.13371,  -1.15048,   0.19149,   0.45912,   1.07201, $
        1.68603,   0.89978,   1.99886,   2.00232,   1.50626,   1.21505,   1.37760,   1.47085,   1.91271,   2.22154, $
        2.33221,   2.45040,   2.70493,   2.19647,   2.05535,   2.31949,   2.75410,   2.36939,   2.80703,   2.66936  ]

A4 = [  -0.02177,    0.47008,    0.72921,    0.23599,    -0.08688,    1.23882,    0.88928,    -0.28712,   -0.32759,   -0.79537, $
        -1.36658,   -0.50148,   -1.48155,   -1.45969,    -1.05472,   -0.67196,   -0.85054,    -0.96618,   -1.31350,   -1.56074, $
        -1.63082,   -1.89979,   -1.94290,   -1.59554,    -1.39972,   -1.59101,   -2.02287,    -1.64367,   -2.05247,   -1.94981  ]

A5 = [   0.04642,    -0.08760,   -0.19775,   -0.05303,    0.06487,    -0.45990,   -0.26462,    0.12298,    0.09850,    0.23982, $
         0.44572,     0.11220,    0.44119,    0.42864,    0.30570,     0.14381,    0.21706,    0.26384,    0.37295,    0.44630, $
         0.45959,     0.59943,    0.55999,    0.47378,    0.38845,     0.44155,    0.59338,    0.46056,    0.60221,    0.57715  ]




;creating abscissa vector
XX = INDGEN(10000) ; ABSCISSA
LAMBDA_TABLE = NM * 10.0 ; RECALCULATING WAVELENGHT TO ANGSTROMS

;spline interpolation for coefficients
SPL_A0 = SPLINE(LAMBDA_TABLE, A0, XX, SPLINE_TENSION)
SPL_A1 = SPLINE(LAMBDA_TABLE, A1, XX, SPLINE_TENSION)
SPL_A2 = SPLINE(LAMBDA_TABLE, A2, XX, SPLINE_TENSION)
SPL_A3 = SPLINE(LAMBDA_TABLE, A3, XX, SPLINE_TENSION)
SPL_A4 = SPLINE(LAMBDA_TABLE, A4, XX, SPLINE_TENSION)
SPL_A5 = SPLINE(LAMBDA_TABLE, A5, XX, SPLINE_TENSION)


;calculation of limb darkening factor using formula for wavelenght eq LAMBDA and angle eq THETA
DARKENING_FACTOR =  SPL_A0(LAMBDA) $
              + SPL_A1(LAMBDA)*COS(THETA)^1. $
              + SPL_A2(LAMBDA)*COS(THETA)^2. $
              + SPL_A3(LAMBDA)*COS(THETA)^3. $
              + SPL_A4(LAMBDA)*COS(THETA)^4. $
              + SPL_A5(LAMBDA)*COS(THETA)^5.


IF KEYWORD_SET(GRAPH) THEN BEGIN



       WINDOW, 0, XSIZE = 1280, YSIZE=600, XPOS=0, YPOS=0, TITLE='COEFFICIENTS'
       LOADCT, 39

       PLOT, LAMBDA_TABLE, A0, XRANGE=[3800, 9000], YSTYLE=1,YRANGE=[-3, 3], PSYM=2, SYMSIZE = 0.5
       OPLOT, XX, SPL_A0, COLOR = 40
       PLOTS, [LAMBDA,LAMBDA],[SPL_A0(LAMBDA)-0.1,SPL_A0(LAMBDA)+0.1], COLOR=40

       OPLOT, LAMBDA_TABLE, A1,PSYM=2, SYMSIZE=0.5
       OPLOT, XX, SPL_A1, COLOR = 80
       PLOTS, [LAMBDA,LAMBDA],[SPL_A1(LAMBDA)-0.1,SPL_A1(LAMBDA)+0.1], COLOR=80

       OPLOT, LAMBDA_TABLE, A2, PSYM=2, SYMSIZE=0.5
       OPLOT, XX, SPL_A2, COLOR = 120
       PLOTS, [LAMBDA,LAMBDA],[SPL_A2(LAMBDA)-0.1,SPL_A2(LAMBDA)+0.1], COLOR=120

       OPLOT, LAMBDA_TABLE, A3, PSYM=2, SYMSIZE=0.5
       OPLOT, XX, SPL_A3, COLOR = 160
       PLOTS, [LAMBDA,LAMBDA],[SPL_A3(LAMBDA)-0.1,SPL_A3(LAMBDA)+0.1], COLOR=160

       OPLOT, LAMBDA_TABLE, A4, PSYM=2, SYMSIZE=0.5
       OPLOT, XX, SPL_A4, COLOR = 200
       PLOTS, [LAMBDA,LAMBDA],[SPL_A4(LAMBDA)-0.1,SPL_A4(LAMBDA)+0.1], COLOR=200

       OPLOT, LAMBDA_TABLE, A5, PSYM=2, SYMSIZE=0.5
       OPLOT, XX, SPL_A5, COLOR = 240
       PLOTS, [LAMBDA,LAMBDA],[SPL_A5(LAMBDA)-0.1,SPL_A5(LAMBDA)+0.1], COLOR=240


       WSHOW, 0



       WINDOW, 2, XSIZE = 1280, YSIZE=300, XPOS=0, YPOS=600, TITLE='LIMB DARKENING'

       TH=FINDGEN(158)/100.

       DARKENING_FACTOR_TH =     SPL_A0(LAMBDA) $
              + SPL_A1(LAMBDA)*COS(TH)^1. $
              + SPL_A2(LAMBDA)*COS(TH)^2. $
              + SPL_A3(LAMBDA)*COS(TH)^3. $
              + SPL_A4(LAMBDA)*COS(TH)^4. $
              + SPL_A5(LAMBDA)*COS(TH)^5.

       PLOT, TH, DARKENING_FACTOR_TH, XSTYLE=1, YSTYLE=1, LINESTYLE=2
       PLOTS, [TH(THETA*100),TH(THETA*100)],[DARKENING_FACTOR_TH(THETA*100)-0.1,DARKENING_FACTOR_TH(THETA*100)+0.1], COLOR=200
       WSHOW, 2



ENDIF

RETURN, DARKENING_FACTOR

END


PRO WHITE, THETA, WAVEL, CA_PROF, GRAPH=GRAPH

; purpose:
; calculating Ca H and Ca K profiles for given position on the solar disc (theta)
; including limb darkening

; usage:
; WHITE, (theta(radians), output_lambda, output_profile [,/graph] )

; outputs:
; array of wavelenghts (from 3914.67 to 3999.91 Angstroems with 0.01 step), array of line profile in respect wavelenghts

; important note:
; limb darkening factors from 3929.47 to 3937.68 Angstroems and from 3965.06 to 3972.79 Angstroems
; were taken from White and Suemoto and were interpolated by spline function
; for other wavelenghts (3915,3954,3980,3991,4000) were calculated using function NECKEL
; you have to be aware that there is some kind of in-continuity in calculated profiles

; profiles are consistent in rage: from 3929.47 to 3937.68 A and from 3965.06 to 3972.79 A
; which are marked by vertical line in the plot when you use /graph keyword

; references:
; White, O.R, Suemoto, Z.   SoPh 3, 523 (1968)

; author:
; Maciej Zapior


; Ca K (3933.680), Ca H (3968.490)

LAMBDA = [  3914.67, 3914.82, 3917.75, 3923.84, 3924.72, 3924.83, 3926.23, 3926.84, 3927.03, 3928.52, $
         3928.79, 3929.47, 3931.40, 3931.58, 3931.81, 3932.98 ,3933.40, 3933.52, 3933.680,3933.86, 3933.96, 3934.380, 3935.16, 3935.53, 3936.12, 3936.39, 3936.93, $
         3937.16, 3937.68, 3938.14, 3939.07, 3939.45, 3939.90, 3940.51, 3942.08, 3947.22, 3949.51, $
         3950.65, 3954.16, 3954.29, 3957.39, 3959.08, 3960.03, 3963.54, 3964.09, 3964.40, 3964.68, $
         3965.06, 3965.31, 3966.25, 3967.02, 3967.65, 3968.490,3969.02, 3969.50, 3969.87, 3970.93, 3971.17, $
         3971.60, 3972.79, 3974.11, 3974.97, 3975.63, 3980.42, 3988.17, 3990.84, 3999.91         ]

PROFILE =  [    86.,  86., 81.,    67.,   64.,  63., 54.,    51.,   48.,  36., $
          34.,    28.8,  17., 16.,    14.,   8.90, 6.17,   6.72, 4.1,    5.97,  6.01,    8.99,  13., 14.,    18.,   20.,  23., $
          24.,    27.4, 30.,    39.,   42.,  46., 51.,    59.,   75.,  80., $
          82.,    83.,  83., 78.,    74.,   70.,  49., 46.,    42.,   40., $
          35.8,   32., 23.,    17.,   13.,  4.25,    10.,   14.,  16., 25.,    27., $
                31.,    42.8,   54.,    61.,    66.,    86.,    94.,    95.,    99.         ]



THETA_TAB  =    [   0.000000,     0.643501,       0.927295,     1.04720,     1.15928,     1.26610,     1.36944,     1.42023 ]

NO_LAMBDA = 11

DARK = FLTARR(NO_LAMBDA,8) ;  (LAMBDA_DARK), (THETA),

DARK(0,*)=  NECKEL([ 0.000000,0.643501,0.927295,1.04720,1.15928,1.26610,1.36944,1.42023 ],3915.)

DARK(1,*)=      [  1.000,      0.860,       0.742,        0.677,         0.624,      0.556,       0.497,        0.454     ]
DARK(2,*)=      [  1.00,         0.89,         0.89,         0.79,         0.75,         0.71,         0.65,         0.62 ]
DARK(3,*)=   [   1.000,         0.874,      0.758,       0.696,        0.634,         0.582,      0.524,       0.486  ]

DARK(4,*)=  NECKEL([ 0.000000,0.643501,0.927295,1.04720,1.15928,1.26610,1.36944,1.42023 ],3954.)

DARK(5,*)=   [   1.000,      0.858,       0.728,        0.660,         0.600,      0.534,       0.468,        0.433 ]
DARK(6,*)=   [   1.00,        0.90,        0.86,        0.81,        0.74,        0.70,        0.64,        0.62  ]
DARK(7,*)=   [   1.000,         0.858,      0.728,       0.652,        0.583,         0.515,      0.448,       0.404  ]

DARK(8,*)=  NECKEL([ 0.000000,0.643501,0.927295,1.04720,1.15928,1.26610,1.36944,1.42023 ],3980)

DARK(9,*)=  NECKEL([ 0.000000,0.643501,0.927295,1.04720,1.15928,1.26610,1.36944,1.42023 ],3991)


DARK(10,*)= NECKEL([ 0.000000,0.643501,0.927295,1.04720,1.15928,1.26610,1.36944,1.42023 ],4000.)

;ADD ALSO LAMBDAS FOR DARKS, IT MUST HAVE NO_LAMBDA ELEMENTS
LAMBDA_DARK =   [3915.,  3929.47, 3933.68,    3937.68,   3954.,    3965.06,   3968.49,  3972.79, 3980.,3991, 4000.]

;CREATING OUTPUT LAMBDAS
XX = FINDGEN(8525)/100. + 3914.67

TH = FINDGEN(1571)/1000.

DARK_THETA = FLTARR(NO_LAMBDA)

FOR I=0,NO_LAMBDA-1 DO BEGIN
    SPL_DARK = SPLINE(THETA_TAB, DARK(I,*), TH)
    DARK_THETA(I) = SPL_DARK(THETA*1000.)
ENDFOR

PROFILE_SPL_DARK = SPLINE(LAMBDA_DARK, DARK_THETA, XX)
PROFILE_SPL      = SPLINE(LAMBDA, PROFILE, XX)

PROFILE_OUT = PROFILE_SPL_DARK*PROFILE_SPL

IF KEYWORD_SET(GRAPH) THEN BEGIN

    WINDOW, 0, XSIZE = 1280, YSIZE=600, XPOS=0, YPOS=0
       LOADCT, 39

    PLOT, LAMBDA, PROFILE, XSTYLE=1, YRANGE=[0,100],YSTYLE=1, PSYM=2, SYMSIZE=0.5
    OPLOT, XX, PROFILE_OUT
    PLOTS, [3929.47,3929.47],[0,100], COLOR=120
    PLOTS, [3937.68,3937.68],[0,100], COLOR=120
    PLOTS, [3965.06,3965.06],[0,100], COLOR=120
    PLOTS, [3972.79,3972.79],[0,100], COLOR=120

ENDIF

WAVEL = XX
CA_PROF = PROFILE_OUT

;stop

END


FUNCTION CONTINUUM, LAMBDA, SPLINE_TENSION, GRAPH=GRAPH, CGS=CGS

; purpose:
; calculate solar spectrum continuum for given wavelenght
;
; usage:
; result = continuum (wavelenght [, spline_tension, /GRAPH, /CGS])

; input:
; wavelenght in angstoems

; output:
; intensity of solar continuum in  egr * sr^(-1) * cm^(-2) * A^(-1) * s^(-1)

; keywords:
; /graph - showing variations of solar continnum and chosen wavelenght
; /cgs - return value in units:    egr * sr^(-1) * cm^(-2) * Hz^(-1) * s^(-1)

; references:
; Allen

; author
; Maciek Zapior


IF N_PARAMS() EQ 0 THEN BEGIN
    STOP, 'USAGE:  result = continuum ( wavelenght [, spline_tension, /GRAPH, /CGS] )'
ENDIF

;table from Allen p.172 column 5:
I = [   0.04, 0.20,   0.30, 0.5,    1.19,  2.15,    2.83,  3.01,    3.20,  $
       3.62, 4.1,    4.4,   4.58, 4.60,   4.59, 4.55,   4.54, 4.48,   $
       4.40, 4.31,   4.08, 3.68,   3.27, 2.88,   2.53, 2.24,   $
       1.97, 1.55,   1.23, 0.99,   0.81, 0.564,  0.403,   0.268,    $
       0.183,    0.081, 0.041,  0.0135,  0.0057   ]

;table from Allen p.172 column 1:
MI =[   0.20, 0.22,   0.24, 0.26,   0.28, 0.30,   0.32, 0.34,   0.36, $
       0.37, 0.38,   0.39, 0.40,   0.41, 0.42,   0.43, 0.44,   0.45, $
       0.46, 0.48,   0.50, 0.55,   0.60, 0.65,   0.70, 0.75,   $
       0.80, 0.9,    1.0,   1.1,  1.2, 1.4,    1.6,   1.8,  $
       2.0,  2.5, 3.0,    4.0,   5.0     ]

LAMBDA = ROUND(LAMBDA)

IF LAMBDA LT 2000 OR LAMBDA GT 50000 THEN BEGIN
    PRINT, 'INPUT VALUE MUST BE BETWEEN 2000 AND 50000'
    STOP
ENDIF

if N_PARAMS() EQ 1 THEN BEGIN
    SPLINE_TENSION = 10.0
    ENDIF
XX = FINDGEN(50001)/10000.

I_SPL = SPLINE(MI, I, XX, SPLINE_TENSION)

IF KEYWORD_SET(GRAPH) THEN BEGIN
    LOADCT,38
    WINDOW, 0, XSIZE=1280, YSIZE=640
    PLOT, MI*10000., I, PSYM=2, XRANGE = [2000,9000], XSTYLE=1
    OPLOT, XX*10000., I_SPL, COLOR=150
    PLOTS, [LAMBDA,LAMBDA],[I_SPL(LAMBDA)-0.5, I_SPL(LAMBDA)+0.5], COLOR=200
    WSHOW, 0

ENDIF

IF KEYWORD_SET(CGS) THEN BEGIN

    VALUE = I_SPL(LAMBDA) * LAMBDA^2./3.
    FACTOR = 10.^6 / 10.^18
    RETURN, VALUE * FACTOR

ENDIF ELSE BEGIN

    RETURN, I_SPL(LAMBDA) * 10.^6

ENDELSE

END


FUNCTION DARKENING, THETA, LAMBDA, SPLINE_TENSION, GRAPH=GRAPH

;   purpose:
; calculate limb darkening of the Sun using Allen's formula
;
;    usage:
; RESULT = DARKENING ( theta(radians), wavelenght(angstroems) [, spline tension, /GRAPH] )
;

; output:
; factor of limb darkening for given wavelenght in respect to disc center

; keywords:
; /graph - showing graphs of: 1. v2 coefficient, 2. u2 coefficient, 3. variation of limb darkening for given theta

; References:
; Allen
;
; author:
; Maciek Zapior

IF (N_PARAMS() EQ 0) OR (N_PARAMS() EQ 1)  THEN BEGIN
    PRINT, 'USAGE:  RESULT = DARKENING ( theta ( radians ),   wavelenght ( angstroems )   [, spline tension, /GRAPH] )'
    STOP
ENDIF

if N_PARAMS() EQ 2 THEN SPLINE_TENSION = 1.0

; table taken from Allen:

MI= [   0.20, 0.22,   0.245,   0.265,    0.28,       $
       0.30, 0.32,   0.35, 0.37,   0.38,        $
       0.40, 0.45,   0.50, 0.55,   0.60, 0.80,   $
       1.0,  1.5, 2.0,    3.0,   5.0,  10.0 ]

U = [   0.12, -1.3,   -0.1,    -0.1,  0.38,       $
        0.74, 0.88,   0.98, 1.03,   0.92,          $
        0.91, 0.99,   0.97,    0.93,  0.88,    0.73,  $
        0.64, 0.57,   0.48, 0.35,   0.22, 0.15    ]

V = [   0.33, 1.6,    0.85,  0.90,    0.57,         $
       0.20, 0.03,   -0.10,    -0.16, -0.05,       $
       -0.05,    -0.17, -0.22,  -0.23,   -0.23,    -0.22,     $
       -0.20,    -0.21, -0.18,  -0.12,   -0.07,    -0.07 ]

;creating abscissa vector
XX = INDGEN(10000) ; ABSCISSA
LAMBDA_TABLE = MI * 10.0^4 ; RECALCULATING WAVELENGHT TO ANGSTROMS

;spline interpolation for U and V coefficients
SPL_V = SPLINE(LAMBDA_TABLE, V, XX, SPLINE_TENSION)
SPL_U = SPLINE(LAMBDA_TABLE, U, XX, SPLINE_TENSION)

;calculation of limb darkening factor using Allen's formula for wavelenght eq LAMBDA and angle eq THETA
DARKENING_FACTOR = 1.0 - SPL_U(LAMBDA) - SPL_V(LAMBDA) + SPL_U(LAMBDA)*COS(THETA) + SPL_V(LAMBDA)*(COS(THETA))^2.0

IF KEYWORD_SET(GRAPH) THEN BEGIN

    WINDOW, 0, XSIZE = 1280, YSIZE=300, XPOS=0, YPOS=0, TITLE='V2 COEFFICIENT'
    SPL_V = SPLINE(LAMBDA_TABLE, V, XX, SPLINE_TENSION)
    PLOT, LAMBDA_TABLE, V, XRANGE=[3800, 9000], YSTYLE=1, YRANGE=[-0.25, 0.25],PSYM=2
    OPLOT, XX, SPL_V, COLOR = 150
    PLOTS, [LAMBDA,LAMBDA],[SPL_V(LAMBDA)-0.1,SPL_V(LAMBDA)+0.1], COLOR=200
    WSHOW, 0


    WINDOW, 1, XSIZE = 1280, YSIZE=300, XPOS=0, YPOS=300, TITLE='U2 COEFFICIENT'
    SPL_U = SPLINE(LAMBDA_TABLE, U, XX, SPLINE_TENSION)
    PLOT, LAMBDA_TABLE, U, XRANGE=[3800, 9000], YSTYLE=1, YRANGE=[0.68, 1.08], PSYM=2
    OPLOT, XX, SPL_U, COLOR = 150
    PLOTS, [LAMBDA,LAMBDA],[SPL_U(LAMBDA)-0.1,SPL_U(LAMBDA)+0.1], COLOR=200
    WSHOW, 1

    WINDOW, 2, XSIZE = 1280, YSIZE=300, XPOS=0, YPOS=600, TITLE='LIMB DARKENING'

    TH=FINDGEN(158)/100.

    DARKENING_FACTOR_TH = 1.0 - SPL_U(LAMBDA) - SPL_V(LAMBDA) + SPL_U(LAMBDA)*COS(TH) + SPL_V(LAMBDA)*(COS(TH))^2.0
    PLOT, TH, DARKENING_FACTOR_TH, XSTYLE=1, YSTYLE=1, LINESTYLE=2
    PLOTS, [TH(THETA*100),TH(THETA*100)],[DARKENING_FACTOR_TH(THETA*100)-0.1,DARKENING_FACTOR_TH(THETA*100)+0.1], COLOR=200
    WSHOW, 2



ENDIF

RETURN, DARKENING_FACTOR

END


PRO PEPIK, thetarad, lambdas, signals, LINE_INDEX, SPLINE_TENSION, GRAPH=GRAPH

; purpose:
; calculate David's profile for Halpha or Hbeta
; or Heinzel's profiles for Ca (3933.680 A) or Ca (8542.120 A)
; for quiet solar photosphere

; usage:
; dawid, thetarad(radians), array for wavelenghts, array for line signal profile, line index [, spline tension, /GRAPH]
; line index: 0=none, 1=Halpha, 2=Ca (3933.680), 3=none, 4=Hbeta, 5=Ca (8542.120)
; spline tension - see IDL help for SPLINE function
; /GRAPH keyword - to draw reference David's profiles with calculated profile

; outputs:
; lambdas - vector of wavelenghts of profile in angstoems
; signals - vector of intensities of profile in percents of solar continuum in center of the disc

; note:
; for thetarad greater than 1.4293 (sin(theta)=0.99) procedure gives values the same as for 1.4293
; because of the limit of David;s profiles
;
; procedure uses function DARKENING.PRO which is used to calculation limb darkening for different wavelenght and theta

; references:
; David, K.-H.
; Heinzel

; author:
; Maciek Zapior

IF N_PARAMS() NE 4 AND N_PARAMS() NE 5 THEN BEGIN
    PRINT, '---------------------------------------------------------'
    PRINT, 'USAGE:'
    PRINT, 'PEPIK, thetarad(radians), array for wavelenghts, array for line signal profile, line index [, spline tension, /GRAPH]'
    PRINT, 'line index: 0=none, 1=Halpha, 2=Ca H (3933.680), 3=He D3, 4=Hbeta, 5=Ca (8542.120), 6=Ca K (3933.680)'
    PRINT, 'spline tension - see IDL help for SPLINE function'
    PRINT, '/GRAPH keyword - to draw reference David profiles with calculated profile'
    PRINT, '---------------------------------------------------------'
    STOP
ENDIF

if N_PARAMS() EQ 4 THEN SPLINE_TENSION = 15.0

;ROUND TO 0.0001
THETARAD = ROUND(THETARAD*10000.)/10000.

;PREVENT TO ERRORS
IF THETARAD GT 1.4293 THEN THETARAD = 1.4293

; d3 lambda: 5875.989, 587.5618, 587.49, 5875,7  ??

LAMBDA_CENTR = [0.0,    6562.808,  3968.492,    5875.618,  4861.342,    8542.144,  3933.680]
NO_LAMBDA    = [0.0,       28,      43,        51,      28,        50,     43     ] ; FOR D3 - NO REVERSABLE

; line index: 0=none, 1=Halpha, 2=Ca (3933.680), 3=D3 He, 4=Hbeta, 5=Ca (8542.120)

CASE LINE_INDEX OF

    0: BEGIN
       PRINT,'---------------------------------------------------'
       PRINT, 'YOU HAVE TO CHOOSE LINE INDEX: 1=Halpha, 2=Ca H (3933.680), 3=He D3, 4=Hbeta, 5=Ca (8542.120)'
       PRINT,'---------------------------------------------------'
       STOP
       END


    1: BEGIN
       DELTA_LAMBDA = [  0.0, 0.1,    0.2,   0.3,  0.4, 0.5,    0.6,   0.7,  0.8, 0.9,    $
                   1.0, 1.5,    2.0,   2.5,  3.0, 3.5,    4.0,   4.5,  5.0, 6.0,    $
                   7.0, 8.0,    9.0,   10.0, 15.0,   20.0, 25.0,   30.0            ]

       END


    2: BEGIN
       DELTA_LAMBDA = [  0.000,   0.025,   0.050,   0.075,   0.100,   0.125,   0.150,   0.175,   0.200,   0.225, $
                     0.250,   0.275,   0.300,   0.325,   0.350,   0.375,   0.400,   0.425,   0.450,   0.475, $
                     0.500,   0.525,   0.550,   0.575,   0.600,   0.800,   1.000,   1.200,   1.400,   1.600, $
                     1.800,   2.000,   2.200,   2.400,   2.600,   2.800,   3.000,   3.200,   3.400,   3.600, $
                     3.800,   4.000,   4.200                                           ]
       END


    3: BEGIN
       DELTA_LAMBDA = FINDGEN(51) * 0.5 - 12.5
       ;DELTA_LAMBDA = [ 0.0,    0.5,   1.0,  1.5, 2.0,    2.5,   3.0,  3.5, 4.0,    4.5,   $
       ;           5.0,    5.5,   6.0,  6.5, 7.0,    7.5,   8.0,  8.5, 9.0,    9.5       ]
       DELTA_LAMBDA = DELTA_LAMBDA - 0.118
       END


    4: BEGIN
       DELTA_LAMBDA = [  0.0, 0.1,    0.2,   0.3,  0.4, 0.5,    0.6,   0.7,  0.8, 0.9,    $
                   1.0, 1.5,    2.0,   2.5,  3.0, 3.5,    4.0,   4.5,  5.0, 6.0,    $
                   7.0, 8.0,    9.0,   10.0, 15.0,   20.0, 25.0,   30.0            ]
       END


    5: BEGIN
       DELTA_LAMBDA = [  0.000,   0.020,    0.040, 0.060,  0.080,   0.100,  0.120,    0.140, 0.160,  0.180,  $
                   0.200,  0.220,  0.240,  0.260,  0.280,  0.300,  0.320,  0.340,  0.360,  0.380,   $
                   0.400,  0.420,  0.440,  0.460,  0.480,  0.520,  0.560,  0.600,  0.640,  0.680,   $
                   0.720,   0.760,  0.800,  0.840,  0.880,  0.920,  0.960,  1.000,  1.100,  1.200,    $
                   1.300,  1.400,  1.500,  1.600,  1.700,  1.800,  1.900,  2.000,  2.200,  2.400    ]
       END

    6: BEGIN
       DELTA_LAMBDA = [  0.000,   0.025,   0.050,   0.075,   0.100,   0.125,   0.150,   0.175,   0.200,   0.225, $
                     0.250,   0.275,   0.300,   0.325,   0.350,   0.375,   0.400,   0.425,   0.450,   0.475, $
                     0.500,   0.525,   0.550,   0.575,   0.600,   0.800,   1.000,   1.200,   1.400,   1.600, $
                     1.800,   2.000,   2.200,   2.400,   2.600,   2.800,   3.000,   3.200,   3.400,   3.600, $
                     3.800,   4.000,   4.200                                           ]
       END

ENDCASE



; profil


THETA_DAVID  = [    0.0,   0.6435,   0.9273,   1.1198,   1.2532,   1.3705,   1.4293]  ; RADIANS


LINE_PROF=fltarr(NO_LAMBDA(LINE_INDEX),7)   ; (POSITION IN WAVELENGHT, POSITION ON THE DISC - THETA)


CASE LINE_INDEX OF
    1: BEGIN

         ;HA PROFILE
         LINE_PROF(*,0)=[ 16.9,   17.5, 18.9,   22.8, 30.2,   40.4, 50.6,   56.8, 61.0,   63.8, $
                    66.0,   72.5, 77.0,   80.3, 82.9,   84.9, 86.6,   88.0,     89.3,  90.9,   $
                    92.2,   93.2, 94.3,   95.0, 97.3,   98.4, 99.0,   99.3 ]

         LINE_PROF(*,1)=[ 17.0,   17.8, 19.6,   23.5, 30.6,   41.0, 51.2,   59.7, 65.0,   67.9, $
                    70.1,   76.7, 80.8,   83.6, 85.8,   87.6, 89.1,   90.3,     91.3,  92.9,   $
                    94.2,   95.1,    95.9,     96.6,  98.2,   98.9,    99.4,     99.5 ]

         LINE_PROF(*,2)=[ 18.1,   18.6, 20.2,   23.9, 30.9,   41.3, 53.3,   62.2, 68.2,   72.0, $
                    74.4,   80.5, 84.8,   87.7, 89.7,   91.0, 92.2,   93.1,     93.9,  95.0,   $
                    96.1,   96.7,    97.3,     97.8,  99.0,   99.5,    99.7,     99.8 ]

         LINE_PROF(*,3)=[ 19.0,   19.4, 20.8,   24.1, 31.2,   41.7, 55.4,   66.4, 72.5,   76.5, $
                    79.0,   85.5, 88.9,   91.1, 92.7,   93.9, 94.9,   95.6,     96.1,  96.9,   $
                    97.6,   98.2,    98.6,     99.0,  99.6,   99.8,    100.0,    100.0 ]

         LINE_PROF(*,4)=[ 20.2,   20.6, 21.9,   25.3, 31.7,   42.8, 56.6,   67.5, 75.1,   80.3, $
                    84.2,   89.9, 92.6,   94.2, 95.6,   96.5, 97.0,   97.4,     97.7,  98.3,   $
                    98.7,   99.0,    99.2,     99.4,  99.8,   99.9,    100.0,    100.0 ]

         LINE_PROF(*,5)=[ 22.6,   23.0, 24.9,   28.0, 34.5,   46.5, 58.7,   72.0, 81.5,   85.5, $
                    88.7,   93.8, 95.5,   96.7, 97.3,   97.9, 98.4,   98.6,     98.9,  99.2,   $
                    99.4,   99.6,    99.7,     99.8,  100.0,  100.0,  100.0,  100.0 ]

         LINE_PROF(*,6)=[ 30.6,   30.9, 31.8,   34.2, 39.0,   48.2, 62.5,   77.1, 87.4,   91.7, $
                    94.1,   96.5, 97.7,   98.3, 98.8,   99.1, 99.3,   99.4,     99.5,  99.7,   $
                    99.8,   99.9,    99.9,     100.0,     100.0,     100.0,     100.0,     100.0 ]
       END
    2: BEGIN

         LINE_PROF(*,0)=[  4.250,  4.482,  4.914,  5.433,  6.025,  6.458,  6.760,  6.864,  6.945,  7.052, $
                     7.209,  7.414,  7.675,  7.931,  8.194,  8.426,  8.710,  8.960,  9.214,  9.455, $
                     9.726,  9.979, 10.196, 10.389, 10.610, 12.266, 14.099, 15.312, 16.371, 17.841, $
                    19.451, 21.111, 22.833, 24.618, 26.425, 28.212, 29.987, 32.075, 34.609, 36.969, $
                    38.926, 40.586, 42.552                                     ]

         LINE_PROF(*,1)=[     4.497,  4.617,  5.031,  5.502,  6.040,  6.700,  6.983,  7.106,  7.184,  7.291, $
                        7.374,  7.621,  7.912,  8.188,  8.442,  8.655,  8.848,  9.131,  9.349,  9.597, $
                        9.799, 10.009, 10.167, 10.376, 10.652, 12.325, 14.196, 15.442, 16.538, 18.058, $
                       19.721, 21.438, 23.222, 25.075, 26.957, 28.825, 30.687, 32.875, 35.527, 38.009, $
                       40.087, 41.864, 43.964                                    ]

         LINE_PROF(*,2)=[     5.306,  5.380,  5.719,  6.190,  6.788,  7.466,  7.959,  8.269,  8.453,  8.571, $
                        8.645,  8.800,  9.013,  9.242,  9.507,  9.765,  9.986, 10.273, 10.517, 10.738, $
                       10.944, 11.161, 11.345, 11.570, 11.890, 13.704, 15.725, 17.040, 18.177, 19.773, $
                    21.513, 23.297, 25.137, 27.037, 28.954, 30.840, 32.708, 34.903, 37.570, 40.038, $
                       42.063, 43.760, 45.776                                    ]

         LINE_PROF(*,3)=[     6.226,  6.312,  6.636,  7.135,  7.809,  8.513,  9.204,  9.874, 10.178, 10.351, $
                       10.417, 10.456, 10.540, 10.746, 11.015, 11.287, 11.543, 11.853, 12.096, 12.362, $
                       12.635, 12.885, 13.117, 13.395, 13.760, 15.787, 18.025, 19.441, 20.633, 22.329, $
                       24.172, 26.046, 27.960, 29.918, 31.871, 33.770, 35.627, 37.813, 40.480, 42.902, $
                       44.822, 46.371, 48.230                                        ]

         LINE_PROF(*,4)=[  7.157,  7.249,  7.522,  7.999,  8.679,  9.492, 10.349, 11.002, 11.552, 11.883, $
                    11.953, 11.969, 12.058, 12.167, 12.309, 12.557, 12.886, 13.200, 13.476, 13.767, $
                    14.013, 14.304, 14.572, 14.932, 15.229, 17.449, 19.908, 21.451, 22.741, 24.585, $
                    26.585, 28.614, 30.682, 32.791, 34.890, 36.927, 38.910, 41.245, 44.101, 46.682, $
                    48.708, 50.325, 52.274                                     ]

         LINE_PROF(*,5)=[  7.996,  8.179,  8.394,  8.793,  9.423, 10.299, 11.237, 12.190, 13.037, 13.577, $
                    13.671, 13.640, 13.639, 13.623, 13.745, 13.883, 14.081, 14.387, 14.663, 14.940, $
                    15.123, 15.486, 15.844, 16.224, 16.682, 19.147, 21.861, 23.581, 25.031, 27.087, $
                    29.322, 31.597, 33.926, 36.309, 38.688, 41.000, 43.259, 45.921, 49.175, 52.129, $
                    54.469, 56.354, 58.622                                     ]

         LINE_PROF(*,6)=[  8.540,  8.658,  8.855,  9.283,  9.886, 10.788, 11.861, 12.744, 13.744, 14.473, $
                    14.720, 14.640, 14.580, 14.529, 14.623, 14.723, 14.794, 15.025, 15.319, 15.552, $
                    15.643, 16.015, 16.433, 16.889, 17.369, 19.981, 22.855, 24.702, 26.276, 28.487, $
                    30.899, 33.367, 35.904, 38.515, 41.133, 43.688, 46.197, 49.156, 52.774, 56.084, $
                    58.739, 60.909, 63.516                                     ]


       END
    3: BEGIN

         LINE_PROF(*,0)=[  98.989, 98.305, 97.393, 98.620, 98.683, 98.534, 99.082, 77.069, 97.275, 92.490, $
                     98.226, 99.004, 98.254, 98.650, 98.815, 99.067, 98.879, 98.803, 98.643, 98.778, $
                     97.602, 98.713, 98.269, 98.803, 98.765, 98.289,                     $
                     98.674, 97.918, 99.017, 98.662, 97.738, 98.860, 98.327, 89.609, 90.006, 98.785, $
                     98.351, 98.565, 98.729, 98.704, 97.985, 98.003, 96.925, 97.358, 96.294, 97.734, $
                     96.513, 97.823, 97.719, 95.197, 95.211                             ]


       END
    4: BEGIN

         ;HB PROFILE

         LINE_PROF(*,0)=[ 13.3,   14.9, 20.3,   30.9, 41.9,   48.2, 51.9,   55.0, 57.9,   60.2, $
                    62.4,   69.5, 73.8,   77.3, 79.6,   81.8, 83.5,   84.8,     86.1,  88.1,   $
                    89.7,   91.1, 92.2,   93.1, 96.0,   97.6, 98.4,   98.9  ]

         LINE_PROF(*,1)=[ 14.1,   15.8, 21.1,   31.9, 45.9,   52.7, 56.9,   60.1, 62.9,   65.1, $
                    67.0,   73.7, 77.7,   81.3, 83.8,   85.9, 87.2,   88.4,     89.3,  91.0,   $
                    92.7,   93.9,    94.8,     95.7,  97.9,   98.7,    99.1,     99.4   ]

         LINE_PROF(*,2)=[ 15.2,   16.7, 22.3,   33.8, 49.2,   57.6, 62.5,   65.7, 68.5,   70.7, $
                    72.7,   79.2, 83.1,   86.0, 88.3,   89.8, 91.0,   92.0,     92.6,  94.0,   $
                    95.0,   96.0,    96.7,     97.3,  98.9,   99.4,    99.6,     99.8   ]

         LINE_PROF(*,3)=[ 16.3,   17.7, 23.7,   36.2, 55.4,   63.3, 68.4,   72.1, 74.9,   77.1, $
                    78.9,   84.8, 88.1,   90.5, 92.3,   93.5, 94.3,   95.0,     95.6,  96.5,   $
                    97.3,   97.8,    98.3,     98.6,  99.4,   99.7,    99.8,     99.9   ]

         LINE_PROF(*,4)=[ 18.0,   19.3, 25.0,   38.0, 60.4,   70.6, 75.8,   79.1, 81.4,   83.1, $
                    84.6,   89.7, 92.3,   94.1, 95.3,   96.1, 96.7,   97.1, 97.5,   98.1,     $
                    98.5,   98.8, 99.0,   99.2, 99.7,   99.8, 99.9,   100.0 ]

         LINE_PROF(*,5)=[ 20.0,   21.1, 26.5,   40.6, 64.5,   76.3, 82.4,   86.0, 88.3,   89.8, $
                    91.1,   94.5, 96.4,   97.4, 97.9,   98.4, 98.6,   98.8, 99.0,   99.2, $
                    99.3,   99.5, 99.7,   99.9, 100.0,  100.0,   100.0,    100.0  ]

         LINE_PROF(*,6)=[ 21.5,   22.6, 28.0,   42.5, 69.5,   81.6, 86.8,   88.7, 91.4,   92.7, $
                    93.7,   96.6, 97.9,   98.6, 98.9,   99.1, 99.3,   99.4, 99.5,   99.6,     $
                    99.7,   99.8, 99.9,   100.0,   100.0,   100.0,   100.0,   100.0    ]
       END
    5: BEGIN
         LINE_PROF(*,0)=[ 17.900, 17.950, 18.300, 18.800, 19.450, 20.500, 21.750, 23.600, 25.400, 27.550, $
                    30.250, 32.300, 34.550, 36.600, 38.250, 39.650, 40.900, 41.900, 42.850, 43.700, $
                    44.400, 45.150, 45.550, 46.200, 46.850, 47.750, 48.700, 49.750, 50.550, 51.250, $
                    51.950, 52.700, 53.600, 54.300, 54.950, 55.650, 56.250, 56.600, 58.144, 60.168, $
                    61.742, 63.065, 64.513, 66.030, 67.481, 68.830, 70.078, 71.261, 73.568, 75.829  ]

         LINE_PROF(*,1)=[ 18.311, 18.519, 18.887, 19.475, 20.190, 21.387 ,22.573, 24.051, 26.046, 28.275, $
                    30.785, 34.044, 36.339, 39.083, 41.048, 42.781, 43.400, 45.049, 45.998, 46.806, $
                    47.641, 48.219, 48.800, 49.485, 50.141, 51.186, 52.311, 53.262, 54.077, 54.813, $
                    55.491, 56.218, 57.041, 58.180, 58.896, 59.571, 60.417, 60.911, 62.353, 63.915, $
                    65.247, 66.483, 67.788, 69.122, 70.403, 71.592, 72.673, 73.688, 75.673, 77.627  ]

         LINE_PROF(*,2)=[ 18.584, 18.812, 19.062, 19.617, 20.285, 21.260, 22.697, 24.501, 26.691, 29.000, $
                    31.320, 34.104, 37.262, 39.938, 42.596, 44.757, 45.645, 47.698, 48.931, 49.894, $
                    50.904, 51.600, 52.369, 53.092, 53.795, 54.971, 56.157, 57.076, 57.953, 58.700, $
                    59.408, 60.105, 60.836, 62.049, 62.759, 63.398, 64.238, 64.795, 66.139, 67.455, $
                    68.644, 69.799, 71.001, 72.215, 73.378, 74.449, 75.407, 76.295, 78.032, 79.747  ]

         LINE_PROF(*,3)=[ 18.704, 18.832, 19.128, 19.652, 20.304, 21.267, 22.648, 24.416, 26.564, 28.806, $
                    31.068, 34.062, 37.065, 40.095, 42.859, 45.448, 47.296, 49.497, 51.179, 52.415, $
                    53.598, 54.609, 55.537, 56.292, 57.065, 58.338, 59.489, 60.431, 61.397, 62.131, $
                    62.902, 63.570, 64.215, 65.216, 65.868, 66.471, 67.115, 67.662, 68.930, 70.204, $
                    71.351, 72.445, 73.592, 74.754, 75.858, 76.864, 77.753, 78.569, 80.160, 81.727  ]

         LINE_PROF(*,4)=[ 18.728, 18.748, 18.968, 19.404, 19.980, 20.717, 21.815, 23.525, 25.276, 27.227, $
                    29.621, 32.607, 35.972, 39.541, 42.609, 45.571, 48.426, 50.631, 52.778, 54.306, $
                    55.638, 57.016, 58.067, 58.846, 59.689, 61.030, 62.105, 63.092, 64.147, 64.859, $
                    65.702, 66.344, 66.918, 67.597, 68.174, 68.744, 69.130, 69.633, 70.845, 72.184, $
                    73.342, 74.393, 75.517, 76.667, 77.750, 78.726, 79.581, 80.362, 81.876, 83.361  ]

         LINE_PROF(*,5)=[ 18.764, 18.671, 18.835, 19.194, 19.706, 20.243, 21.095, 22.753, 24.161, 25.858, $
                    28.365, 31.351, 35.040, 39.020, 42.407, 45.675, 49.473, 51.691, 54.252, 56.039, $
                    57.506, 59.205, 60.370, 61.170, 62.077, 63.480, 64.490, 65.517, 66.650, 67.344, $
                    68.250, 68.868, 69.378, 69.921, 70.388, 70.853, 71.365, 71.724, 72.617, 74.006, $
                    75.171, 76.182, 77.283, 78.419, 79.480, 80.427, 81.250, 81.996, 83.437, 84.843  ]

         LINE_PROF(*,6)=[ 18.783, 18.635, 18.772, 19.092, 19.572, 20.011, 20.741, 22.373, 23.612, 25.183, $
                    27.744, 30.734, 34.586, 38.772, 42.312, 45.745, 50.012, 52.241, 55.008, 56.926, $
                    58.461, 60.321, 61.542, 62.355, 63.292, 64.729, 65.706, 66.752, 67.925, 68.610, $
                    69.548, 70.154, 70.632, 71.140, 71.570, 72.000, 72.460, 72.762, 73.529, 74.941, $
                    76.108, 77.098, 78.187, 79.315, 80.364, 81.296, 82.101, 82.829, 84.231, 85.597  ]
       END
    6: BEGIN
            LINE_PROF(*,0)=[     4.100,  4.201,  4.607,  5.147,  5.640,  6.034,  6.305,  6.297,  6.200,  6.121, $
                           6.074,  6.080,  6.170,  6.270,  6.405,  6.649,  6.825,  7.041,  7.225,  7.430, $
                           7.585,  7.822,  7.979,  8.146,  8.305,  9.572, 10.431, 11.295, 12.307, 13.120, $
                          13.730, 14.923, 16.355, 17.668, 19.066, 20.408, 21.632, 22.739, 23.751, 24.867, $
                          26.133, 27.409, 28.610                                        ]

         LINE_PROF(*,1)=[     4.760,  4.835,  5.091,  5.523,  5.998,  6.456,  6.748,  6.869,  6.822,  6.619, $
                     6.519,  6.565,  6.647,  6.834,  7.031,  7.297,  7.547,  7.827,  8.055,  8.243, $
                     8.471,  8.765,  8.911,  9.079,  9.254, 10.623, 11.534, 12.440, 13.500, 14.338, $
                    14.945, 16.179, 17.661, 19.005, 20.427, 21.779, 22.991, 24.070, 25.040, 26.109, $
                    27.326, 28.543, 29.674                                     ]

         LINE_PROF(*,2)=[     5.786,  5.861,  6.093,  6.497,  6.990,  7.513,  7.887,  8.148,  8.246,  8.163, $
                     8.074,  8.044,  8.044,  8.126,  8.276,  8.478,  8.709,  8.949,  9.188,  9.405, $
                     9.659,  9.967, 10.154, 10.368, 10.551, 12.080, 13.077, 14.059, 15.212, 16.110, $
                    16.741, 18.064, 19.658, 21.089, 22.597, 24.014, 25.269, 26.368, 27.339, 28.409, $
                    29.633, 30.849, 31.959                                     ]

         LINE_PROF(*,3)=[     6.813,  6.949,  7.244,  7.648,  8.099,  8.704,  9.278,  9.678, 10.022, 10.106, $
                    10.031,  9.939,  9.835,  9.812,  9.788,  9.845, 10.005, 10.170, 10.380, 10.632, $
                    10.874, 11.189, 11.383, 11.639, 11.908, 13.615, 14.720, 15.805, 17.079, 18.062, $
                    18.744, 20.197, 21.949, 23.512, 25.155, 26.692, 28.045, 29.219, 30.248, 31.383, $
                    32.681, 33.966, 35.130                                     ]

         LINE_PROF(*,4)=[     7.819,  7.899,  8.122,  8.550,  9.149,  9.873, 10.539, 11.187, 11.596, 11.828, $
                    11.877, 11.754, 11.599, 11.419, 11.285, 11.176, 11.260, 11.361, 11.560, 11.811, $
                    11.997, 12.326, 12.631, 12.879, 13.108, 14.997, 16.227, 17.436, 18.856, 19.954, $
                    20.722, 22.341, 24.297, 26.050, 27.895, 29.624, 31.151, 32.479, 33.649, 34.940, $
                    36.419, 37.886, 39.217                                     ]

         LINE_PROF(*,5)=[      8.733,  8.780,  9.045,  9.466, 10.121, 10.948, 11.729, 12.619, 13.371, 13.918, $
                    14.060, 13.983, 13.781, 13.500, 13.187, 12.983, 12.951, 13.075, 13.214, 13.368, $
                    13.570, 13.946, 14.134, 14.433, 14.616, 16.756, 18.165, 19.559, 21.194, 22.473, $
                    23.384, 25.262, 27.531, 29.583, 31.749, 33.792, 35.614, 37.215, 38.643, 40.219, $
                    42.023, 43.820, 45.468                                     ]

         LINE_PROF(*,6)=[     9.200,  9.251,  9.545,  9.997, 10.706, 11.583, 12.536, 13.550, 14.521, 15.213, $
                         15.515, 15.489, 15.267, 14.959, 14.595, 14.349, 14.250, 14.310, 14.318, 14.380, $
                         14.529, 14.949, 15.155, 15.427, 15.621, 17.934, 19.470, 20.996, 22.787, 24.197, $
                         25.215, 27.281, 29.777, 32.048, 34.450, 36.725, 38.766, 40.575, 42.198, 43.994, $
                         46.045, 48.094, 49.987                                       ]

       END

ENDCASE



;creating abscissa vector
XX = FINDGEN(14294) / 10000.; ABSCISSA

IF LINE_INDEX NE 3 THEN BEGIN

    ;spline interpolation for U and V coefficients
    SPL_LINE = FLTARR(14294,NO_LAMBDA(LINE_INDEX))

    FOR T=0,NO_LAMBDA(LINE_INDEX)-1 DO BEGIN
       SPL_LINE(*,T) = SPLINE(THETA_DAVID, LINE_PROF(T,*), XX, SPLINE_TENSION)
    ENDFOR

    lambdas=fltarr(NO_LAMBDA(LINE_INDEX)*2-1)
    signals=fltarr(NO_LAMBDA(LINE_INDEX)*2-1)

    THETA_INDEX = ROUND(THETARAD*10000)

    signals(NO_LAMBDA(LINE_INDEX)-1:2*NO_LAMBDA(LINE_INDEX)-2) = SPL_LINE(THETA_INDEX,*)
    lambdas(NO_LAMBDA(LINE_INDEX)-1:2*NO_LAMBDA(LINE_INDEX)-2) = DELTA_LAMBDA(*) + LAMBDA_CENTR(LINE_INDEX)

    signals(0:NO_LAMBDA(LINE_INDEX)-2) = reverse(SPL_LINE(THETA_INDEX,1:NO_LAMBDA(LINE_INDEX)-1),2)
    lambdas(0:NO_LAMBDA(LINE_INDEX)-2) = -1.0*reverse(DELTA_LAMBDA(1:NO_LAMBDA(LINE_INDEX)-1))+LAMBDA_CENTR(LINE_INDEX)

ENDIF ELSE BEGIN
    SIGNALS = LINE_PROF(*,0)
    FOR I=1,6 DO LINE_PROF(*,I) = LINE_PROF(*,0)
    LAMBDAS = DELTA_LAMBDA + LAMBDA_CENTR(LINE_INDEX)
    PRINT, 'D3'
    ;STOP
ENDELSE

;USING FUNCTION TO CALCULATE LIMB DARKENING
LIMB_DARK = DARKENING(THETARAD,ROUND (LAMBDA_CENTR(LINE_INDEX)))

signals = signals * LIMB_DARK


IF KEYWORD_SET(GRAPH) THEN BEGIN
    LOADCT,39
    WINDOW,0, XSIZE=1280, YSIZE=640
    PLOT, LAMBDA_CENTR(LINE_INDEX)+DELTA_LAMBDA, LINE_PROF(*,0), XRANGE = [LAMBDA_CENTR(LINE_INDEX)-10,LAMBDA_CENTR(LINE_INDEX)+10], XSTYLE=1

    FOR T=1,6 DO BEGIN
       OPLOT, LAMBDA_CENTR(LINE_INDEX)+ DELTA_LAMBDA, LINE_PROF(*,T)*DARKENING(THETA_DAVID(T) , ROUND (LAMBDA_CENTR(LINE_INDEX))), COLOR=T*255./7.+ 30.
    ENDFOR

    OPLOT, LAMBDAS, SIGNALS, LINESTYLE=2

ENDIF

end


PRO SNAKE
START:
window,1,xpos=750,ypos=0,xsize=300,ysize=300,title='SNAKE PLISSKEN'
XYOUTS, 150, 270, 'SNAKE PLISSKEN ', CHARTHICK=1.5, CHARSIZE=1.7, COLOR=255, ALIGNMENT=0.5,/device
XYOUTS, 150, 245, 'The first ever game for IDL ', CHARTHICK=1.5, CHARSIZE=1.7, COLOR=255, ALIGNMENT=0.5,/device
XYOUTS, 150, 185, 'authors: ', CHARTHICK=1, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.5,/device
XYOUTS, 150, 145, 'Maciej Zapior ', CHARTHICK=1, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.5,/device
XYOUTS, 150, 125, 'Astronomocal Institute of Wroclaw University', CHARTHICK=1, CHARSIZE=0.85, COLOR=255, ALIGNMENT=0.5,/device
XYOUTS, 150, 80, 'Piotr Podgorski ', CHARTHICK=1, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.5,/device
XYOUTS, 150, 60, 'Space Reaserch Centre Polish Academy of Science ', CHARTHICK=1, CHARSIZE=0.85, COLOR=255, ALIGNMENT=0.5,/device
MENU,['PLAY','USAGE'],KK,TITLE='DO YOU READY TO PLAY?',BUTSIZE=350
wdelete,1
IF KK EQ 2 THEN BEGIN
    ERASE
    XYOUTS, 100, 250, 'W - UP ', CHARTHICK=1.5, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.0,/device
    XYOUTS, 100, 200, 'S - DOWN ', CHARTHICK=1.5, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.0,/device
    XYOUTS, 100, 150, 'A - LEFT ', CHARTHICK=1.5, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.0,/device
    XYOUTS, 100, 100, 'D - RIGHT ', CHARTHICK=1.5, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.0,/device
    XYOUTS, 100,  50, 'Q - QUIT ', CHARTHICK=1.5, CHARSIZE=1.5, COLOR=255, ALIGNMENT=0.0,/device
ENDIF
MENU,['MASTA','SOFT MASTA','INTERMEDIATE','EXPERIENCED','BEGINER'],KK,TITLE='CHOOSE SKILL LEVEL:',BUTSIZE=250
case kk of
    1: rozmiar = 100
    2: rozmiar = 200
    3: rozmiar = 300
    4: rozmiar = 400
    5: rozmiar = 500
else:
endcase
czekaj = 0.010
erase
i = 0
C = ''
D = ''
score = 0
dlugosc = 100
x1 = intarr(8000)
y1 = intarr(8000)
L = 0
H = 0
Idx = 10
Idx2 = 9
Idx3 = 12
Idx4 = 11
loadct,39
window,0,xsize=750,ysize=750,title='SNAKE PLISSKEN IDL FIRST EVER GAME',xpos=0
plot, [0,0], [0,0], xrange=[-0.5*rozmiar,0.5*rozmiar], yrange=[-0.5*rozmiar,0.5*rozmiar], /isotropic, $
    title='SNAKE PLISSKEN IDL FIRST EVER GAME',xstyle=1,ystyle=1
rana = RANDOMU(Seed, 2)
rana = rana * (rozmiar - 20) - 0.5 * (rozmiar - 10)
rannax = ROUND(rana(0))
rannay = ROUND(rana(1))
oplot, [rannax],[rannay], symsize=1.0, color = 255, psym = 6, thick=1.0
WHILE 1 DO BEGIN
I = I + 1
CONTROL = GET_KBRD(0)
wait, CZEKAJ
print, 'c przed = ', c
c = STRLOWCASE(c)
IF CONTROL EQ "a" THEN C = CONTROL
IF CONTROL EQ "s" THEN C = CONTROL
IF CONTROL EQ "d" THEN C = CONTROL
IF CONTROL EQ "w" THEN C = CONTROL
IF CONTROL EQ "q" THEN C = CONTROL
IF (C NE "") AND (H EQ 0 ) THEN H = 1
print,'CONTROL=',CONTROL
print,'C      =',c
CASE C OF
    "q":  goto,KONIEC
    "a":   BEGIN
          x1(Idx) = x1(Idx2) - 1
          y1(Idx) = y1(Idx2)
         END
    "s":   BEGIN
          y1(Idx) = y1(Idx2) - 1
          x1(Idx) = x1(Idx2)
         END

    "w":   BEGIN
          y1(Idx) = y1(Idx2) + 1
          x1(Idx) = x1(Idx2)
         END
    "d":   BEGIN
          x1(Idx) = x1(Idx2) + 1
          y1(Idx) = y1(Idx2)
          END
    ELSE:
ENDCASE
    IF x1(Idx) LE -0.5*rozmiar THEN goto,KONIEC
    IF x1(Idx) GE 0.5*rozmiar THEN goto,KONIEC
    IF y1(Idx) LE -0.5*rozmiar THEN goto,KONIEC
    IF y1(Idx) GE 0.5*rozmiar THEN goto,KONIEC
    L = Idx + 2
    L = L mod dlugosc
    for k = 0, (dlugosc - 5) do begin
    IF (x1(L) EQ x1(Idx)) AND (y1(L) EQ y1(Idx)) AND H GT dlugosc THEN goto,KONIEC
    L = L + 1
    L = L mod dlugosc
    endfor
    if (x1(Idx) ge rannax-2 AND x1(Idx) le rannax+2) and (y1(Idx) ge rannay-2 AND y1(Idx) le rannay+2)  then begin
       PRINT,'KAMYK'
       oplot, [rannax],[rannay], symsize=1.0, color = 0, psym = 6, thick=1.0
       rana = RANDOMU(Seed, 2)
       rana = rana * (rozmiar - 20) - 0.5 * (rozmiar - 10)
       rannax = ROUND(rana(0))
       rannay = ROUND(rana(1))
       dlugosc = dlugosc + 50 ;wysypuje sie
       oplot, [rannax],[rannay], symsize=1.0, color = 255, psym = 6, thick=1.0
       score = score + 1
    ENDIF
    IF H GT 0 THEN H = H + 1
    IF H GT (dlugosc + 10) THEN H = (dlugosc + 10)
    oplot, [x1(Idx2),x1(Idx)],[y1(Idx2),y1(Idx)] , symsize=0.5, color = 255, psym = 6, thick=2.5
    oplot, [x1(Idx4),x1(Idx3)],[y1(Idx4),y1(Idx3)] , symsize=0.5, color = 0, psym = 6, thick=2.5
    oplot, [rannax],[rannay], symsize=1.0, color = 255, psym = 6, thick=1.0
    Idx = Idx + 1
    Idx = Idx mod dlugosc
    Idx2 = Idx2 + 1
    Idx2 = Idx2 mod dlugosc
    Idx3 = Idx3 + 1
    Idx3 = Idx3 mod dlugosc
    Idx4 = Idx4 + 1
    Idx4 = Idx4 mod dlugosc
ENDWHILE
KONIEC:
if c eq "q" then begin
    XYOUTS, 0, -10, 'QUIT THE GAME', CHARTHICK=3, CHARSIZE=3.5, COLOR=225, ALIGNMENT=0.5,/DATA
    GOTO, KONIEC_PO_Q
endif
PRINT, 'GAME OVER'
XYOUTS, 400, 375, 'GAME OVER ', CHARTHICK=3, CHARSIZE=3.5, COLOR=225, ALIGNMENT=0.5,/DEVICE
XYOUTS, 405, 260, 'SCORE: ', CHARTHICK=2, CHARSIZE=2.5, COLOR=50, ALIGNMENT=0.5,/DEVICE
XYOUTS, 400, 200, strcompress(string(SCORE*100),/remove_all), CHARTHICK=2, CHARSIZE=2.5, COLOR=50, ALIGNMENT=0.5,/DEVICE
KONIEC_PO_Q:
MENU,['YES','NO'],KK,TITLE='DO YOU WANT START AGAIN?',BUTSIZE=350
IF KK EQ 1 THEN GOTO, START
END










;-----------------------------------------------------------------------------------------------------------------------------------------------------------
;---------   SPECTRAL ANALISYS  MODULE   -----------------------------------------------------------------------------------------------------------------
;-----------------------------------------------------------------------------------------------------------------------------------------------------------

;reading file with spectrum lines

CLOSE,/ALL
ON_IOERROR, c223
A1=' '
OPENR,1,'lines.txt'
FOR I=0.,999. DO READF,1,A1
c223: no_ref_lines=I
CLOSE,1
ON_IOERROR,NULL


openr,1,'lines.txt'


;NO_REF_LINES=15

line = {name:'',wave:0.0,inten:0.0,atm_color:0}
lines = replicate(line,NO_REF_LINES)

line_by_line = strarr(NO_REF_LINES)

for iii=0,no_ref_lines-1 do begin
    A=''
    readf,1,FORMAT = '(A16,F8.3,F7.3,I4)', A,B,C,D
    LINES(III).NAME=A
    LINES(III).WAVE=B
    LINES(III).INTEN=C
    LINES(III).ATM_COLOR=D
endfor

close,1

TOTAL_BEGINNING:

menu,['SJ','SJ HA','SJ HA CA','SJ HA CA D3 ','SJ HA CA D3 HB','SJ HA CA D3 HB IR','OTHER SELECTION'],NO_CAM,TITLE='SELECT CAMERAS',BUTSIZE=240

cam_index = strarr(6)


LITT = strarr(6)
LITT(0) = 'sj'
LITT(1) = 'ha'
LITT(2) = 'ca'
LITT(3) = 'd3'
LITT(4) = 'hb'
LITT(5) = 'ir'

IF NO_CAM EQ 7 THEN BEGIN
    TEMP_INDEX = 0

    repeat begin

       MENU,$
       ['SJ',$
        'HA',$
        'CA',$
        'D3',$
        'HB',$
        'IR',$
        '.:: OK ::.'],$
       LETR,TITLE='SELECT CAMERAS (one click for one camera)',BUTSIZE=635,YOFFSET=50

       IF LETR LT 7 THEN CAM_INDEX(TEMP_INDEX) = litt(letr-1)

       IF LETR NE 7 THEN TEMP_INDEX = TEMP_INDEX + 1

    ENDREP UNTIL LETR EQ 7
    NO_CAM = TEMP_INDEX

ENDIF ELSE BEGIN

    cam_index(0) = 'sj'
    cam_index(1) = 'ha'
    cam_index(2) = 'ca'
    cam_index(3) = 'd3'
    cam_index(4) = 'hb'
    cam_index(5) = 'ir'

ENDELSE



;SCAT_spectra=fltarr(1280,no_cam)

;MENU,['MEAN','MEDIAN'],calculation_method,TITLE='CALCULATE METHOD',BUTSIZE=160,xoffset=0,yoffset=0

calculation_method=1

REVERSE_CONTROL = BYTARR(NO_CAM)
REF_TIME_CONTROL = BYTARR(NO_CAM)

DISPERSION_COEF = DBLARR(2,NO_CAM)

WAVELENGTH_DATA = FLTARR(1280,NO_CAM)
;FACTOR_DATA = DBLARR(NO_CAM) ; CALIBRATION INTENSITY FACTOR FOR EACH LINE
FACTOR_DATA_VECTOR = FLTARR(1280,NO_CAM)
FACTOR_DATA_COEF = FLTARR(3,NO_CAM)
FACTOR_DATA_COEF(0,*) = 1.0  ; NO INTENSITY CALIBRATION, OTHER COEFFICIENTS EQ 0.0
FACTOR_CGS = DBLARR(NO_CAM)  ; FACTOR TO TRANSFER % CONTINUUM TO CGS UNITS
FACTOR_EXP = DBLARR(NO_CAM)  ; FACTOR TO CORRECT EXPOSITION TIMES RATIO

REF_EXP_TIME = DBLARR(NO_CAM)

FOR I=0,NO_CAM-1 DO FACTOR_DATA_VECTOR(*,I) = 1.
FOR I=0,NO_CAM-1 DO FACTOR_CGS(I)  = 1.
FOR I=0,NO_CAM-1 DO FACTOR_EXP(I)  = 1.
FOR I=1,NO_CAM-1 DO WAVELENGTH_DATA(*,I)=DINDGEN(1280)

SUN_X_CENTR = -1. ; IT MEANS THAT X,Y COORDINATES OF DISC CENTER ARE NOT CALCULATED
SUN_Y_CENTR = -1.

REV=INTARR(NO_CAM) ; INFORMATION WHICH SPECTRUM WAS REVERSED 0=NO, 1=YES
REV=0

MIN_WL = FLTARR(6)
MAX_WL = FLTARR(6)

	;for Ca K line
	MIN_WL = [0.,6553.,3919.,5863.,4848.,8529.5]
	MAX_WL = [0.,6573.,3946.,5888.,4875.,8555.5]
	LAMBDA_CENTR = [0.0, 6562.808, 3933.682, 5875.618, 4861.342, 8542.144  ]


SUN_RADIUS = 3159.*2. ;;;  <<<---------------------------------------------------  check it  -------------------------------------------

menu,['MANUALLY','AUTO'],LETR_MAN,TITLE='FINDING SPECTRAL IMAGES',BUTSIZE=600





;---------   START OF THE LOOP IMAGE BY IMAGE  --------------------

START_FILE:

CGS_DONE = BYTARR(NO_CAM)

CALIBRATION_COMPLETE = 0
SUN_CENTRE_CONTROL = 0
Y_TITLE = STRARR(NO_CAM)



CASE LETR_MAN OF

1: BEGIN  ;manual finding

    ;menu,['SJ','SJ HA','3','4','5','6'],LETR,TITLE='NUMBER OF CAMERAS',BUTSIZE=600

    data_name = strarr(NO_CAM)

    data_name(0) = dialog_pickfile(title='CHOOSE SJ IMAGE', /file,  /MULTIPLE_FILES, get_path = dir, FILTER = ['*sj*_R.fit'])

;   cam_index = strarr(6)
;   cam_index(0) = 'sj'
;   cam_index(1) = 'ha'
;   cam_index(2) = 'ca'
;   cam_index(3) = 'd3'
;   cam_index(4) = 'hb'
;   cam_index(5) = 'ir'

    FOR I=1,NO_CAM-1 DO BEGIN

       data_name(I)=dialog_pickfile(title='CHOOSE'+cam_index(i)+ 'SPECTRAL IMAGE', /file,  path = dir, FILTER = ['*'+cam_index(i)+'*_R.fit'])

    ENDFOR

    ;no_cam = size(file_name,/n_elements)

    print,data_name

;-------------------
    name_length=strlen(data_name)
    slash_pos=strpos(data_name,'\', /REVERSE_SEARCH ) ; because whole path
    rev_slash_pos = name_length - slash_pos
    write_name=strmid(data_name,slash_pos+1, rev_slash_pos)
    write_len =strlen(write_name)
    write_name=strmid(write_name,0, write_len)
    write_name=strcompress(write_name)
    print,'write_name ',write_name(0)
    data_time = strmid(write_name(0),9,10)
    print,'data_time ',data_time
    ;print,data_name(0)

;----------------
    EXP_TIME = FLTARR(NO_CAM)
    FOR I=0,NO_CAM-1 DO BEGIN
       ;data_name(I) = findfile(dir+'*'+data_time(0)+'_'+cam_index(i)+'*_R.fit')
       N_pos=strpos(data_name(I),'N', /REVERSE_SEARCH )
       R_pos=strpos(data_name(I),'R' ,/REVERSE_SEARCH )
       EXP_LEN = R_POS-N_POS
       TEMP = strmid(data_name(I),N_POS+2,EXP_LEN)
       EXP_TIME(I) = FLOAT(TEMP)
    ENDFOR


    END



2: BEGIN  ;auto finding

    data_name = strarr(NO_CAM)
    data_name(0) = dialog_pickfile(title='CHOOSE SJ IMAGE', /file,  /MULTIPLE_FILES, get_path = dir, FILTER = ['*sj*_R.fit'])
;   cam_index = strarr(6)
;   cam_index(0) = 'sj'
;   cam_index(1) = 'ha'
;   cam_index(2) = 'ca'
;   cam_index(3) = 'd3'
;   cam_index(4) = 'hb'
;   cam_index(5) = 'ir'

;-------------------

;    20080515_063239.000_ca_F_30.fit

    name_length=strlen(data_name)
    slash_pos=strpos(data_name,'\', /REVERSE_SEARCH ) ; because whole path
    rev_slash_pos = name_length - slash_pos
    write_name=strmid(data_name,slash_pos+1, rev_slash_pos)
    write_len =strlen(write_name)
    write_name=strmid(write_name,0, write_len)
    write_name=strcompress(write_name)
    print,'write_name ',write_name(0)
    data_time = strmid(write_name(0),9,10)
    print,'data_time ',data_time
    ;print,data_name(0)

;----------------
    EXP_TIME = FLTARR(NO_CAM)
    FOR I=0,NO_CAM-1 DO BEGIN
       data_name(I) = findfile(dir+'*'+data_time(0)+'_'+cam_index(i)+'*_R.fit')
       N_pos=strpos(data_name(I),'N', /REVERSE_SEARCH )
       R_pos=strpos(data_name(I),'R' ,/REVERSE_SEARCH )
       EXP_LEN = R_POS-N_POS
       TEMP = strmid(data_name(I),N_POS+2,EXP_LEN)
       EXP_TIME(I) = FLOAT(TEMP)
    ENDFOR
;----------------

    ;print,write_name
    ;stop

    END
ENDCASE

;0000000000000000000000000000000000000000000000000000000000000000000000000000
;0000000000000000000000000000000000000000000000000000000000000000000000000000
;0000000000000000000000000000000000000000000000000000000000000000000000000000



;CONTOL IF FILE WITH COEFFICIENTS EXIST
FILE_CONTROL = FINDFILE('calibration_coefficients.txt')

;READING FILE WITH COEFFICIENTS
IF FILE_CONTROL(0) EQ 'calibration_coefficients.txt' THEN BEGIN

    PRINT, 'READING calibration_coefficients.txt FILE '

    CLOSE,/ALL
    ON_IOERROR, ERROR_1
    A1=' '
    OPENR,1,'calibration_coefficients.txt'
    FOR I=0.,999. DO READF,1,A1
    ERROR_1: NO_LINES=I
    CLOSE,1
    ON_IOERROR,NULL

    IF NO_LINES NE (NO_CAM*2+NO_CAM*3+NO_CAM) THEN BEGIN
       PRINT, 'CHECK NUMBER OF CAMERAS ANALYSED AND NUMBER OF COEFFICIENT IN calibration_coefficients.txt FILE'
       PRINT, 'NUMBER OF COEFFICIENT SHOULD BE EQUAL 6 TIMES NUMBER OF CAMERAS'
       PRINT, 'YOU CAN DELETE FILE calibration_coefficients.txt AND START ANALYSIS FROM THE BEGINNING AND MAKE CALIBRATION ONE MORE TIME'
       MENU,['OK'],LINE_CORR,TITLE='YOU HAVE TO ANALYSE ALL AND ONLY CAMERAS WHICH YOU CALIBRATE BEFORE',BUTSIZE=800,xoffset=0,yoffset=-15
       MENU,['OK'],LINE_CORR,TITLE='YOU CAN DELETE FILE calibration_coefficients.txt ',BUTSIZE=800,xoffset=0,yoffset=-15
       MENU,['OK'],LINE_CORR,TITLE='IN THIS CASE START ANALYSIS FROM THE BEGINNING AND MAKE CALIBRATION ONE MORE TIME',BUTSIZE=800,xoffset=0,yoffset=-15



       GOTO,TOTAL_BEGINNING
    ENDIF


    OPENR,1,'calibration_coefficients.txt'
    ;readf,1,FORMAT = '(A16,F8.3,F7.3,I4)', A,B,C,D
    ;READF,1,SLIT_WIDTH
    WHOLE_FILE = DBLARR(NO_CAM*2+NO_CAM*3+NO_CAM)
    READF,1,WHOLE_FILE
    CLOSE,1

    K=0
    FOR I=0,NO_CAM-1 DO BEGIN
       FOR J=0,1 DO BEGIN
         DISPERSION_COEF(J,I) = WHOLE_FILE(K)
         K=K+1
       ENDFOR
    ENDFOR

    FOR I=0,NO_CAM-1 DO BEGIN
       FOR H=0,2 DO BEGIN
         FACTOR_DATA_COEF(H,I) = WHOLE_FILE(K)
         K=K+1
       ENDFOR
    ENDFOR

    FOR I=0,NO_CAM-1 DO BEGIN
       REF_EXP_TIME(I) = WHOLE_FILE(K)
       K=K+1
    ENDFOR

    FOR I=1,NO_CAM-1 DO BEGIN
       WAVELENGTH_DATA(*,I)=FINDGEN(1280.)*DISPERSION_COEF(1,I) + DISPERSION_COEF(0,I)
       ;FACTOR_EXP(I) = REF_EXP_TIME(I) / EXP_TIME(I)
       Y_TITLE(I) = 'PERCENT OF SOLAR CONTINUUM'
    ENDFOR

    CALIBRATION_COMPLETE = 1
    REF_TIME_CONTROL(*) = 1

ENDIF



;CORRECTION OF EXPOSURE TIMES RATIO:
;FOR I=1, NO_CAM-1 DO BEGIN
;   IF REF_TIME_CONTROL(I) EQ 1 THEN BEGIN
;     FACTOR_EXP(I) = REF_EXP_TIME(I) / EXP_TIME(I)
;   ENDIF
;ENDFOR
;STOP

data = fltarr(1280,512,no_cam)
hair_pos_up_coef_data = fltarr(2,no_cam)
hair_pos_down_coef_data = fltarr(2,no_cam)


; reading SJ image
data(*,*,0)=readfits(data_name(0),header_SJ)

hair_pos_up_coef_data  (0,0) =  float(strcompress(strmid(header_SJ( 7),9,21),/remove_all))
hair_pos_down_coef_data(0,0) =  float(strcompress(strmid(header_SJ( 8),9,21),/remove_all))

IF N_ELEMENTS(header_SJ) GT 11 THEN BEGIN
    SUN_X_CENTR =           float(strcompress(strmid(header_SJ( 9),9,21),/remove_all))
    SUN_Y_CENTR =           float(strcompress(strmid(header_SJ(10),9,21),/remove_all))
    IF N_ELEMENTS(header_SJ) GT 13 THEN BEGIN
       X_SLIT_POSITION =        float(strcompress(strmid(header_SJ(11),9,21),/remove_all))
    ENDIF
ENDIF

;DETERMINING SIZE OF HEADER OF SPECTRAL IMAGES ON THE EXAMPLE OF HA IMAGE

data(*,*,1)=readfits(data_name(1),header)

;CREATING ARRAY FOR HEADERS:
;HEADER_SPECTRA = STRARR(N_ELEMENTS(HEADER),NO_CAM)
HEADER_SPECTRA = STRARR(20,NO_CAM)

;reading spectral images
for i=1,no_cam-1 do begin
    data(*,*,i)=readfits(data_name(i),header)
    hair_pos_up_coef_data(0,i) =          float(strcompress(strmid(header( 7),9,21),/remove_all))
    hair_pos_down_coef_data(0,i) =        float(strcompress(strmid(header( 8),9,21),/remove_all))

    IF N_ELEMENTS(header) GT 11 THEN BEGIN

       DISPERSION_COEF (0,I) =    float(strcompress(strmid(header( 9),9,21),/remove_all))
       DISPERSION_COEF (1,I) =    float(strcompress(strmid(header(10),9,21),/remove_all))

       WAVELENGTH_DATA(*,I)=FINDGEN(1280.)*DISPERSION_COEF(1,I) + DISPERSION_COEF(0,I)
       ;FACTOR_EXP(I) = REF_EXP_TIME(I) / EXP_TIME(I)

       IF N_ELEMENTS(header) GT 13 THEN BEGIN

         ;FAC_DEGREE =           float(strcompress(strmid(header(11),9,21),/remove_all))
         X_COORD = FINDGEN(1280)
         FACTOR_DATA_VECTOR(*,I) = 0.

         FOR H=0,2 DO BEGIN

          FACTOR_DATA_TEMP =      float(strcompress(strmid(header(11+H),9,21),/remove_all))
          FACTOR_DATA_VECTOR(*,I) = FACTOR_DATA_VECTOR(*,I) + X_COORD^H*FACTOR_DATA_TEMP
          PRINT, FACTOR_DATA_TEMP
          FACTOR_DATA_COEF(H,I) =  FACTOR_DATA_TEMP

         ENDFOR

         Y_TITLE(I) = 'PERCENT OF SOLAR CONTINUUM'

       ENDIF
    ENDIF

    FOR H=0,N_ELEMENTS(header)-1 DO BEGIN
       HEADER_SPECTRA(H,I) = HEADER(H)
    ENDFOR
endfor



;line_pos_coef(0,0) is position of the slit on the SJ image
openr,1,dir+'..\resources\'+'coefficients_line.txt'
readf,1,line_pos_coef
close,1

X_SLIT_POSITION = line_pos_coef(0,0)

;000000000000000000000000
;hair_pos_up_coef_data(0,1) =   hair_pos_up_coef_data(0,1) +4.
;hair_pos_down_coef_data(0,1) = hair_pos_down_coef_data(0,1)+4.
;00000000000000000000000000

;---  PARAMETERS OF PROGRAM ---

;SIZE OF IMAGES
xsize=1280
ysize=512

;ENLARGEMENT OF OUTPUT IMAGES TO MOVE INPUT IMAGES IN THE INTERIOR OF THEM
;number of pixels added on edges
xedge=0
yedge=500

;reading images:

;stop

for i=0,no_cam-1 do wdelete,i

image = dblarr(xsize+2*xedge,ysize+2*yedge,no_cam)

IMAGE(xedge:xsize+xedge-1,yedge:ysize+yedge-1,*) = data

;calculating current width
hair_width = dblarr(no_cam)
hair_width     = abs(hair_pos_up_coef_data(0,*) - hair_pos_down_coef_data(0,*))

hair_ref_width = abs(hair_pos_up_coef_data(0,0) - hair_pos_down_coef_data(0,0))

;scale of transformation
scale = dblarr(no_cam)
scale = double(hair_ref_width) / double(hair_width)

;scaling image
scaled_image =      dblarr(xsize+2*xedge,ysize+2*yedge,no_cam)
;scaled_image_10 =  dblarr(xsize+2*xedge,(ysize+2*yedge)*10,no_cam)


;calculation of the limits of spectra on sj image
down_limit = fltarr(no_cam)
up_limit  =  fltarr(no_cam)

;calculation of the limit of spectra in coordinates of SJ image
;it is further shown as dashed horizontal lines with caption "end if HA spectrum"
for i=1,no_cam-1 do begin
    down_limit(i) = hair_pos_down_coef_data(0,0) - scale(i)*hair_pos_down_coef_data(0,I)
    up_limit(i) = hair_pos_up_coef_data(0,0) - scale(i)*(hair_pos_up_coef_data(0,I)-511.)
endfor


;new width of whole image
new_ysize = dblarr(no_cam)
new_ysize = round(scale(*) * ysize)
;new_ysize_10 = round(scale(*) * ysize * 10.)

for i=0,no_cam-1 do begin
;   IF SCALE GE 1 THEN BEGIN
       scaled_image (xedge:xsize+xedge-1,  yedge:new_ysize(i)+yedge-1,i) = $
       congrid(image(xedge:xsize+xedge-1,  yedge:ysize+yedge-1       ,i),xsize,new_ysize(i),1,/interp,/minus_one)
endfor

;stop

new_hair1 = dblarr(no_cam)
new_hair2 = dblarr(no_cam)

translate_1 = dblarr(no_cam)
translate_2 = dblarr(no_cam)
translate   = dblarr(no_cam)

for i=0,no_cam-1 do begin

    ;new_hair1(i) = float(hair_pos_down_coef_data(0,i)-yedge)*scale(i)+yedge
    ;new_hair2(i) = float(hair_pos_up_coef_data(0,i)-yedge)*scale(i)+yedge

    translate_1(i) = hair_pos_down_coef_data(0,i)*scale(i) - hair_pos_down_coef_data(0,0)
    translate_2(i) = hair_pos_up_coef_data(0,i)*scale(i)   - hair_pos_up_coef_data(0,0)

    ;to safety:
    translate(i) = (translate_1(i) + translate_2(i)) * 0.5
    translate(i) = round(translate(i))

    print,'SCALE: ',scale(I), '  TRANSLATE: ',translate(i)

endfor


;moving image
trans_image = fltarr(xsize+2*xedge,ysize+2*yedge,no_cam)

;for i=1,no_cam-2 do begin

;   trans_image (xedge : xsize+xedge-1 , yedge-translate(i) : new_ysize(i)+yedge-translate(i)-1,i) = $
;   scaled_image(xedge : xsize+xedge-1 , yedge              : new_ysize(i)+yedge-1,i)

;endfor

for i=0,no_cam-1 do begin

    trans_image (*,*,I) = shift(scaled_image(*,*,i),0,-translate(i))
    PRINT,I,TRANSLATE(I)

endfor

 ;trans_image(*,*,0) = scaled_image(*,*,0)

for i=0,no_cam-1 do wdelete,i

;cutting edges in respect to sj limit of image:

trans_image_cut = fltarr(1280,512,no_cam)
trans_image_cut(*,*,*) = trans_image(*,yedge:yedge+511,*)

;REVERSE SPECTRAL IMAGES !!!
;REMEMBER IF THEY ARE REVERSE OR NOT

;IMPORTANT NOTE:
;THERE IS NO REVERSE OF IR SPECTRUM !!!

FOR I=1,NO_CAM-1 DO BEGIN
    IF I NE 5 THEN BEGIN
       TRANS_IMAGE_CUT(*,*,I) = REVERSE(TRANS_IMAGE_CUT(*,*,I))
    ENDIF
ENDFOR
;wait, 0.5

for j=0,no_cam-1 do begin

       window,j,ypos=0*150,xpos=j*75,xsize=320*4,ysize=128*4,title=cam_index(j)+' .:: AFTER RE-SCALING AND SHIFTING ::. '
       ;tvscl,congrid(trans_image(*,*,j),1380/1.,612/1.)
       tvscl,trans_image_cut(*,*,j)
       wait,0.4
endfor

WSHOW,0

wait,0.5






FOR I=0,NO_CAM-1 DO FACTOR_CGS(I) = 1.

;IF CALIBRATION_COMPLETE EQ 1 THEN Y_TITLE = 'PERCENT OF SOLAR CONTINUUM'

;XLOADCT,BLOCK=1

;stop


UP_LEVEL =4095.




TEMP_FOR_INTEGRATION = TRANS_IMAGE_CUT



FIRST_AGAIN = 0


;--------------   START LOOP BY REGION OF INTEREST   ------------



START:

IF FIRST_AGAIN NE 0 THEN BEGIN

    window,0,ypos=0*150,xpos=0*75,xsize=320*4,ysize=128*4,title=' .:: SELECT REGION OF INTEREST ::. ( '+data_name(0)+' )'
    tvscl,trans_image_cut(*,*,0) < UP_LEVEL

    PLOTS, [0,1279],[UP_LIMIT(1),UP_LIMIT(1)],/DEVICE,COLOR=255
    PLOTS, [0,1279],[UP_LIMIT(1),UP_LIMIT(1)],/DEVICE,COLOR=0,LINESTYLE=2

    PLOTS, [0,1279],[DOWN_LIMIT(1),DOWN_LIMIT(1)],/DEVICE,COLOR=255
    PLOTS, [0,1279],[DOWN_LIMIT(1),DOWN_LIMIT(1)],/DEVICE,COLOR=0,linestyle=2

    XYOUTS,1104,UP_LIMIT(1)+3,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 0 ; shadow
    XYOUTS,1102,UP_LIMIT(1)+5,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 255 ; white letters

    XYOUTS,1084,DOWN_LIMIT(1)-13,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 0
    XYOUTS,1082,DOWN_LIMIT(1)-11,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 255

    XYOUTS,5,UP_LIMIT(1)+3,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 0
    XYOUTS,3,UP_LIMIT(1)+5,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 255

    XYOUTS,5,DOWN_LIMIT(1)-13,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 0
    XYOUTS,3,DOWN_LIMIT(1)-11,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 255

    CURSOR, x_1, Y, /device, /down

    y_pos_1 = y
    print,y
    PLOTS,[0,1279],[y,y],/device,color=255

    CURSOR, x_2, Y, /device, /down

    y_pos_2 = y
    print,y
    PLOTS,[0,1279],[y,y],/device,color=255

    y_down=y_pos_1<y_pos_2
    y_up  =y_pos_2>y_pos_1

    X_POSITION = (X_1 + X_2) / 2.

     PRINT,'UP',Y_UP
     PRINT,'DOWN',Y_DOWN

    SLIT_WIDTH = Y_UP - Y_DOWN

    Y_POSITION = (Y_UP + Y_DOWN) / 2.

AFTER_SCAT_LIGHT_REMOVED:

    spectra=fltarr(1280,no_cam)

    CASE calculation_method OF

    1: BEGIN

       for j=1,no_cam-1 do begin
         for i=0,1279 do begin
          sum_col = total(trans_image_cut(i,y_down:y_up,j))
          spectra(i,j) = sum_col / (y_up-y_down+1)
         endfor
       endfor

       END


    2: BEGIN

       for j=1,no_cam-1 do begin
         for i=0,1279 do begin
          spectra(i,j) = median(trans_image_cut(i,y_down:y_up,j))
         endfor
       endfor

       END

    ENDCASE

ENDIF



;-----------------------------------------------------------------------------------------------------------------------------------------------------------
;-----------------------------------------------------------------------------------------------------------------------------------------------------------
;-----------------------------------------------------------------------------------------------------------------------------------------------------------


REPEAT BEGIN

IF FIRST_AGAIN NE 0 THEN BEGIN

    ;---  DISPLAYING ANALIZED SPECTRA:
    !P.MULTI = [0,0,no_cam-1,0,0]
    window,no_cam,ypos=80,xpos=80,xsize=320*4,ysize=128*7,title=' SPECTRUM / SLIT WIDTH ='$
    +STRCOMPRESS(STRING(SLIT_WIDTH))+'  /  X= '+STRCOMPRESS(STRING(FIX(X_POSITION)))+'  /  Y= '+STRCOMPRESS(STRING(FIX(Y_POSITION)))
    FOR J=1, NO_CAM-1 DO BEGIN
       PLOT,WAVELENGTH_DATA(*,J),FACTOR_EXP(J)*FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), $
       XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
    ENDFOR
    !P.MULTI = 0

ENDIF



MENU,['EXTRACT SPECTRUM','ANOTHER FILE','QUIT','XLOADCT',$
       'CALCULATION SCATERED LIGTH','CALCULATE DISC CENTER',$
       'REF. PHOTOSPH. PROFILE','LIVE PROFILE','WAVELENGTH CORRECTION',$
       'SUBTRACTION SCATERED LIGTH','DOPPLER VELOCITY','SAVE COEFFICIENTS',$
       'CREDITS','CGS AND INTEGRAL','SHOW OTHER IMAGES','SAVE SPECTRA'],LETR,TITLE='AGAIN?',BUTSIZE=250,xoffset=-50,yoffset=538

    IF (LETR EQ 1) THEN BEGIN
       FIRST_AGAIN = 1
       GOTO,START
    ENDIF
    IF (LETR EQ 2) THEN GOTO,START_FILE
    IF (LETR EQ 4) THEN BEGIN
       XLOADCT,BLOCK=1
    ENDIF
    IF (LETR EQ 5) THEN BEGIN

         START_SCAT:

          REPEAT BEGIN
              MENU,['1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20'],NO_SCAT,TITLE='NUMBER OF POINTS',BUTSIZE=300,xoffset=0,yoffset=-15

              window,0,ypos=0*150,xpos=0*75,xsize=320*4,ysize=128*4,title=' .:: SELECT REGION OF SCATERED LIGTH ::. '
              tvscl,trans_image_cut(*,*,0) < UP_LEVEL

              PLOTS, [0,1279],[UP_LIMIT(1),UP_LIMIT(1)],/DEVICE,COLOR=255,LINESTYLE=2
              PLOTS, [0,1279],[DOWN_LIMIT(1),DOWN_LIMIT(1)],/DEVICE,COLOR=255,LINESTYLE=2
              XYOUTS,1104,UP_LIMIT(1)+3,   'TOP LIMIT OF HA SPECTRUM',/DEVICE
              XYOUTS,1084,DOWN_LIMIT(1)-11,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE
              XYOUTS,5,UP_LIMIT(1)+3,   'TOP LIMIT OF HA SPECTRUM',/DEVICE
              XYOUTS,5,DOWN_LIMIT(1)-11,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE

              SCAT_HALF =  1
              SCAT_WIDTH = 3

              Y_POS_SCAT = DBLARR(NO_SCAT)
              SCAT_spectra = DBLARR(1280,NO_CAM,NO_SCAT)

              FOR K=0,NO_SCAT-1 DO BEGIN

                 CURSOR, x, Y, /device, /down
                 y_pos_SCAT(K) = y
                 print,y

                 PLOTS,[0,1279],[y-SCAT_HALF,y-SCAT_HALF],/device,color=255 ;PLOTTING HORIZONTAL LINES
                 PLOTS,[0,1279],[y+SCAT_HALF,y+SCAT_HALF],/device,color=255

                 PRINT,'POSITION',Y

                 CASE calculation_method OF
                 1: BEGIN
                   for j=1,no_cam-1 do begin
                    for i=0,1279 do begin
                        sum_col = total(trans_image_cut(i,y-SCAT_HALF:y+SCAT_HALF,j))
                        SCAT_spectra(i,j,K) = sum_col / 3. ;(PIXEL ROW,CAM,SCAT POINT)
                    endfor
                   endfor
                   END

                 2: BEGIN  ; TO POKI CO NIE JEST ZMIENIONE
                   for j=1,no_cam-1 do begin
                    for i=0,1279 do begin
                        SCAT_spectra(i,j) = median(trans_image_cut(i,y_down:y_up,j))
                    endfor
                   endfor
                   END
                 ENDCASE

              ENDFOR

              !P.MULTI = [0,0,NO_SCAT,0,0]
              FOR CAM=1, NO_CAM-1 DO BEGIN
                 window,CAM,ypos=0*150+10*CAM,xpos=0*75+10*CAM,xsize=320*4,ysize=128*7,title=' .::'+strupcase(CAM_INDEX(CAM))+' MATRIX OF SCATERED LIGTH SPECTRUM ::. '
                 ;window,J,ypos=80+10*J,xpos=80+10*J,xsize=320*4,ysize=128*7,title=strupcase(CAM_INDEX(J))+' SPECTRUM OF SCATTERED LIGHT'

                 FOR K=0,NO_SCAT-1 DO BEGIN
                   PLOT,SCAT_SPECTRA(*,CAM,K), XSTYLE = 1,TITLE=' SCATTERED LIGHT SPECTRUM AT Y POSITION='+STRING(y_pos_SCAT(K)),charsize=1.5
                 ENDFOR
              ENDFOR
              !P.MULTI = 0

              MENU,['AGAIN','OK','CANCEL'],PYTKO,TITLE='SCATTERED LIGHT',BUTSIZE=200,xoffset=0,yoffset=0
              IF PYTKO EQ 1 THEN GOTO,START_SCAT
              IF PYTKO EQ 3 THEN BREAK



              FOR NN=1,NO_CAM-1 DO BEGIN
                 OPENW,4,'scat_light_profiles_'+strlowcase(CAM_INDEX(NN))+'.txt',/APPEND,WIDTH=70000
                 OPENW,5,'scat_light_position_'+strlowcase(CAM_INDEX(NN))+'.txt',/APPEND
                   FOR JJ=0,NO_SCAT-1 DO BEGIN
                    PRINTF,4,SCAT_SPECTRA(*,NN,JJ)
                    DISTANCE_R = SQRT( (SUN_X_CENTR-X_SLIT_POSITION)^2. + (SUN_Y_CENTR-y_pos_SCAT(JJ)*2.)^2. ) / SUN_RADIUS ; binning
                    PRINTF,5,X_SLIT_POSITION,y_pos_SCAT(JJ),DISTANCE_R
                   ENDFOR
                 CLOSE,4
                 CLOSE,5
              ENDFOR

          ENDREP UNTIL PYTKO EQ 2


    ENDIF

    IF (LETR EQ 6) THEN BEGIN
       MENU,['DISC CENTRE','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20'],NO_POINTS,TITLE='CHOOSE NUMBER OF POINTS',BUTSIZE=250,xoffset=0,yoffset=-15

       IF NO_POINTS EQ 1 THEN BEGIN ; ASSUME THAT OBSERVATIONS WERE MADE AT DISC CENTRE
         SUN_X_CENTR = 640.
         SUN_Y_CENTR = 256.
         SUN_CENTRE_CONTROL = 1
       ENDIF ELSE BEGIN




         window,0,ypos=0*150,xpos=0*75,xsize=320*4,ysize=128*4,title=' .:: MARK SOLAR EDGE ::. '
         tvscl,trans_image_cut(*,*,0) < UP_LEVEL


         NO_POINTS = NO_POINTS + 3
         EDGE_POS_X = DBLARR(NO_POINTS)
         EDGE_POS_Y = DBLARR(NO_POINTS)
         I=0
         window,0,ypos=0*150,xpos=0*75,xsize=320*4,ysize=128*4,title=' .:: MARK SOLAR EDGE ::. '
         tvscl,trans_image_cut(*,*,0) < UP_LEVEL
         REPEAT BEGIN
          CURSOR, X, Y, /device, /down
          PLOTS,[x,x],[Y-2,Y+3],/device,color=255 & PLOTS,[x-2,x+3],[Y,Y],/device,color=255
          EDGE_POS_X(I) = X
          EDGE_POS_Y(I) = Y
          I=I+1
         ENDREP UNTIL I EQ NO_POINTS



         ;here it is usage of  MPFITELLIPSE function
         ;it fits ellipse to vector of data
         ;in our case from clicking on the solar edge on the image
         ;here we use it with /circular key it means than circle is fitted instead of ellipse
         ;two first parameters are equal
         ;also parinfo keyword allows us to fix estimated solar radius (in pixels) during fitting
         ;this radius is varying from approx. 3107 pxl (aphelion) to 3212 pxl (peryhelion) with mean value 3212 pxl
         ; tutorial of function:   http://cow.physics.wisc.edu/~craigm/idl/mpfittut.html

         ;first approximation of fitting
         sample_EDGE_COEFF = MPFITELLIPSE(edge_pos_x,edge_pos_y*2,/circular)

         ;creating sructure of parameters
         pi = replicate({fixed:0, limited:[0,0], limits:[0.D,0.D]},5)

         ;fixing first parameter
         PI(0).FIXED=1

         SOLAR_RADIUS_CALC = SUN_RADIUS ; calculated solar radius in primary focus +/- 50 pxl, so it is too small (max 1.6%) so sol_rad = 3159.
                             ; given above valuse
         CHROM_PHOT_RATIO = 1.0060 ;WAS : 1.0028 -> arbitrary value


         MENU,['PHOTOSPHERE','CHROMOSPHERE'],REF_EDGE,TITLE='CHOOSE EDGE',BUTSIZE=300,xoffset=0,yoffset=-15
         IF REF_EDGE EQ 1 THEN BEGIN
          SOL_RAD =  SOLAR_RADIUS_CALC ; IN PIXELS
         ENDIF

         IF REF_EDGE EQ 2 THEN BEGIN
          SOL_RAD =  SOLAR_RADIUS_CALC *  CHROM_PHOT_RATIO; IN PIXELS INCLUDING CHROMOSPHERE HEIGHT
         ENDIF

         ;fitting with fixed value of radius and starting parameters equal to firts approximation:
         EDGE_COEFF = MPFITELLIPSE(edge_pos_x,edge_pos_y*2,[          SOL_RAD,$
                                                  SOL_RAD,$
                                        sample_EDGE_COEFF(2),$
                                        sample_EDGE_COEFF(3),$
                                        sample_EDGE_COEFF(4)],/circular,PARINFO=PI) ;Y*2 BECAUSE OF BINNING in y direction

         SUN_RADIUS=EDGE_COEFF(0)

         IF REF_EDGE EQ 2 THEN BEGIN
          SUN_RADIUS=EDGE_COEFF(0) / CHROM_PHOT_RATIO
         ENDIF


         PLOTS, CIRCLE_BIN(EDGE_COEFF(2), EDGE_COEFF(3), SUN_RADIUS), /Device, COLOR=255
         PLOTS, CIRCLE_BIN(EDGE_COEFF(2), EDGE_COEFF(3), SUN_RADIUS), /Device, COLOR=0, linestyle=2
         SUN_X_CENTR=EDGE_COEFF(2)
         SUN_Y_CENTR=EDGE_COEFF(3)

         ;PRINT,'PREV  RADIUS', PREV_RADIUS
         PRINT,'SOLAR RADIUS', SUN_RADIUS
       ENDELSE
    ENDIF



    IF (LETR EQ 7) THEN BEGIN ;REF PHOTOSPH SPECTRUM

       window,0,ypos=0*150,xpos=0*75,xsize=320*4,ysize=128*4,title=' .:: MARK REFERENCE AREA ON THE SLIT::. '
       tvscl,trans_image_cut(*,*,0) < UP_LEVEL

       PLOTS, [0,1279],[UP_LIMIT(1),UP_LIMIT(1)],/DEVICE,COLOR=255
       PLOTS, [0,1279],[UP_LIMIT(1),UP_LIMIT(1)],/DEVICE,COLOR=0,LINESTYLE=2

       PLOTS, [0,1279],[DOWN_LIMIT(1),DOWN_LIMIT(1)],/DEVICE,COLOR=255
       PLOTS, [0,1279],[DOWN_LIMIT(1),DOWN_LIMIT(1)],/DEVICE,COLOR=0,linestyle=2

       XYOUTS,1104,UP_LIMIT(1)+3,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 0 ; shadow
       XYOUTS,1102,UP_LIMIT(1)+5,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 255 ; white letters

       XYOUTS,1084,DOWN_LIMIT(1)-13,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 0
       XYOUTS,1082,DOWN_LIMIT(1)-11,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 255

       XYOUTS,5,UP_LIMIT(1)+3,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 0
       XYOUTS,3,UP_LIMIT(1)+5,   'TOP LIMIT OF HA SPECTRUM',/DEVICE, color = 255

       XYOUTS,5,DOWN_LIMIT(1)-13,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 0
       XYOUTS,3,DOWN_LIMIT(1)-11,'BOTTOM LIMIT OF HA SPECTRUM',/DEVICE, color = 255

       IF SUN_CENTRE_CONTROL NE 1 THEN BEGIN
         PLOTS, CIRCLE_BIN(SUN_X_CENTR, SUN_Y_CENTR, SUN_RADIUS), /Device, color=0
         ;POSITION OF XYOUTS IN RELATION TO SOLAR EDGE (PLOTTING XYOUTS IN THE CORNER ABOVE SOLAR EDGE)


         IF SUN_X_CENTR LT 0 THEN BEGIN
          POS_X = 1070
         ENDIF ELSE BEGIN
          POS_X = 2
         ENDELSE

         IF SUN_Y_CENTR LT 0 THEN BEGIN
          POS_Y = 476
         ENDIF ELSE BEGIN
          POS_Y = 2
         ENDELSE

         XYOUTS,POS_X,POS_Y+20,'  SOLID LINE - SOLAR EDGE',/DEVICE
         XYOUTS,POS_X,POS_Y,'DASHED LINE - 99% SOLAR RADII',/DEVICE

         PLOTS, CIRCLE_BIN(SUN_X_CENTR, SUN_Y_CENTR, SUN_RADIUS*0.99),LINESTYLE=2, /Device, color=0
         ; 0.99 SOLAR RADII (DAWID GIVES PROFILES UP TO THETA EQ 0.99)

         CURSOR, X, Y, /device, /down
         PLOTS,[0,1279],[y,y],/device,color=255 ;PLOTTING HORIZONTAL LINES
         Y_1 = Y

         CURSOR, X, Y, /device, /down
         PLOTS,[0,1279],[y,y],/device,color=255 ;PLOTTING HORIZONTAL LINES
         Y_2 = Y

         Y_U = Y_1 > Y_2
         Y_D = Y_1 < Y_2
         Y = (Y_D + Y_U) / 2.
         X = X_SLIT_POSITION

         DISTANCE = SQRT((X-SUN_X_CENTR)^2.+(Y*2.-SUN_Y_CENTR)^2.)/SUN_RADIUS ;Y*2 BECAUSE OF BINNING AND ;DIVIDED BY RADIUS
         PRINT,'DISTANCE:',DISTANCE

       ENDIF ELSE BEGIN

         CURSOR, X, Y, /device, /down
         PLOTS,[0,1279],[y,y],/device,color=255 ;PLOTTING HORIZONTAL LINES
         Y_1 = Y

         CURSOR, X, Y, /device, /down
         PLOTS,[0,1279],[y,y],/device,color=255 ;PLOTTING HORIZONTAL LINES
         Y_2 = Y

         Y_U = Y_1 > Y_2
         Y_D = Y_1 < Y_2

         DISTANCE = 0.
       ENDELSE

       REF_WIDTH = Y_U - Y_D + 1

       WAIT,1.0

       REF_spectra=DBLARR(1280,no_cam)

       CASE calculation_method OF
       1: BEGIN
         for j=1,no_cam-1 do begin
          for i=0,1279 do begin
              sum_col = total(trans_image_cut(i,y_D:y_U,j))
              REF_spectra(i,j) = sum_col / (REF_WIDTH)
          endfor
         endfor
         END
       2: BEGIN
         for j=1,no_cam-1 do begin
          for i=0,1279 do begin
              REF_spectra(i,j) = median(trans_image_cut(i,y_D:Y_U,j))
          endfor
         endfor
         END
       ENDCASE

       ;       MENU,['HA','CA','D3','HB','IR'],LINE_REF,TITLE='CHOOSE CAMERA?',BUTSIZE=300,xoffset=0,yoffset=-15
       MENU,[STRUPCASE(CAM_INDEX[1:NO_CAM-1])],LINE_REF,TITLE='CHOOSE CAMERA?',BUTSIZE=300,xoffset=0,yoffset=-15

       CASE CAM_INDEX(LINE_REF) OF
         'ha': LINE_REF_LIEGE = 1
         'ca': LINE_REF_LIEGE = 2
         'd3': LINE_REF_LIEGE = 3
         'hb': LINE_REF_LIEGE = 4
         'ir': LINE_REF_LIEGE = 5
       ENDCASE



       FACTOR_DATA_VECTOR(*,LINE_REF) = 1.

       MED_REF = MEDIAN(REF_spectra(*,LINE_REF))
       FACTOR = MED_REF/100. ; FIRST ESTIMATION OF CALIBRATION FACTOR

;oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo


       IF LINE_REF_LIEGE NE 2 THEN BEGIN
          PEPIK, ASIN(DISTANCE),WAVELENGTH_REF,INTENSITY_REF,LINE_REF_LIEGE
       ENDIF ELSE BEGIN
          WHITE, ASIN(DISTANCE),WAVELENGTH_REF,INTENSITY_REF
          PRINT, 'WHITE USED'
          ;STOP
       ENDELSE

       DARKENING_FACTOR = NECKEL(ASIN(DISTANCE), LAMBDA_CENTR(LINE_REF_LIEGE))



       CLOSE,/ALL
       ON_IOERROR, A323
       A1=' '
       OPENR,1,'liege.txt'
       FOR I=0.,999999. DO READF,1,A1
       A323: no_lines=I
       CLOSE,1
       ON_IOERROR,NULL

       WHOLE_SPECTRUM = DBLARR(2,NO_LINES)
       OPENR,1,'liege.txt'
       READF,1,WHOLE_SPECTRUM
       CLOSE,1

       WAVELENGTH = WHOLE_SPECTRUM(0,*)
       INTENSITY  = WHOLE_SPECTRUM(1,*)

       LOADCT, 39


       WAVE_CORR:

       REPEAT BEGIN

         window,NO_CAM,ypos=80,xpos=80,xsize=320*4,ysize=128*7,title=' REFERENCE SPECTRUM OF PHOTOSPHERE IN '+ STRING(DISTANCE) +'   SOLAR RADII'
         PLOT,WAVELENGTH_DATA(*,LINE_REF),REF_spectra(*,LINE_REF),  XSTYLE=1, TITLE=strupcase(CAM_INDEX(LINE_REF))+' REF SPECTRUM ',charsize=2.0

         IF LINE_REF NE 3 THEN OPLOT, WAVELENGTH_REF, INTENSITY_REF*FACTOR, LINESTYLE=0, COLOR = 150         ; oplot PEPIK profile
         OPLOT, WAVELENGTH,DARKENING_FACTOR*FACTOR*INTENSITY, LINESTYLE=0, COLOR=220  ; oplot LIEGE profile

         MENU,['OK','+10','+1.0','+0.1','+0.01','-0.01','-0.1','-1.0','-10'],CORR,TITLE='WAVELENGTH CORRECTION',BUTSIZE=300,xoffset=0,yoffset=-15
         ;VECTOR OF VALUES;
         CORR=CORR-1
         FAC = [0.,10.,1.0,0.1,0.01,-0.01,-0.1,-1.0,-10.]
         WAVELENGTH_REF = WAVELENGTH_REF + FAC(CORR)
         print, WAVELENGTH_REF(0)

       ENDREP UNTIL CORR EQ 0

       REPEAT BEGIN

         window,NO_CAM,ypos=80,xpos=80,xsize=320*4,ysize=128*7,title=' REFERENCE SPECTRUM OF PHOTOSPHERE IN '+ STRING(DISTANCE) +'   SOLAR RADII'
         PLOT,WAVELENGTH_DATA(*,LINE_REF),REF_spectra(*,LINE_REF),  XSTYLE=1, TITLE=strupcase(CAM_INDEX(LINE_REF))+' REF SPECTRUM ',charsize=2.0


         IF LINE_REF NE 3 THEN OPLOT, WAVELENGTH_REF, INTENSITY_REF*FACTOR, LINESTYLE=0, COLOR = 150  ; oplot PEPIK/DAVID profile

         OPLOT, WAVELENGTH,DARKENING_FACTOR*FACTOR*INTENSITY, LINESTYLE=0,  COLOR=220     ; oplot LIEGE  profile

         ;oplot continuum level = 100. arbitrary units:
         PLOTS, [WAVELENGTH_DATA(0,LINE_REF),WAVELENGTH_DATA(1279,LINE_REF)], [DARKENING_FACTOR*FACTOR*100.,DARKENING_FACTOR*FACTOR*100.], COLOR = 170, LINESTYLE = 2
         ;STOP
         ;OPLOT, WAVELENGTH_REF, (INTENSITY_REF*MED_REF)/100.

         MENU,['OK','+10','+1.0','+0.1','+0.01','-0.01','-0.1','-1.0','-10','WAVELENGTH CORRETION AGAIN'],CORR,TITLE='INTENSITY CORRECTION',BUTSIZE=300,xoffset=0,yoffset=-15
         ;VECTOR OF VALUES;
         CORR=CORR-1
         IF CORR EQ 9 THEN GOTO, WAVE_CORR
         FAC = [0.,10.,1.0,0.1,0.01,-0.01,-0.1,-1.0,-10.]
         FACTOR = FACTOR + FAC(CORR)
         print, FACTOR

       ENDREP UNTIL CORR EQ 0

       MENU,['1','2','3','4','5','6'],NO_AREA,TITLE='HOW MANY AREAS ?',BUTSIZE=300,xoffset=0,yoffset=-15



       IF LINE_REF NE 3 THEN BEGIN   ; CHOOSING PROFILES TO CALIBRATION; EXCEPT D3 LINE; THEN ONLY LIEGE POSSIBLE
         MENU,['LIEGE (RED)','DAVID (GREEN)'],PROFILE_METHOD,TITLE='WHICH PROFILE ?',BUTSIZE=300,xoffset=0,yoffset=-15
       ENDIF ELSE BEGIN
         PROFILE_METHOD = 1
       ENDELSE


       ;ADU_INTENSITY = FLTARR(NO_AREA)
       FACTOR_CALCUL = FLTARR(NO_AREA)
       DATA_C = FLTARR(NO_AREA)

       FOR K=0,NO_AREA-1 DO BEGIN

         CURSOR, X, Y, /DATA, /down
         PLOTS,[X,X],[0,Y],/DATA,color=255 ;PLOTTING VERTICAL LINES
         X_1 = X

         CURSOR, X, Y, /DATA, /down
         PLOTS,[X,X],[0,Y],/DATA,color=255 ;PLOTTING VERTICAL LINES
         X_2 = X

         X_R = X_1 > X_2
         X_L = X_1 < X_2


         ; LIMITS OF AREAS:
         DATA_L = WHERE(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_L) EQ MIN(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_L) ) )
         DATA_R = WHERE(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_R) EQ MIN(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_R) ) )
         DATA_C(K) = (DATA_R + DATA_L) / 2.


         ; LIMITS OF AREAS:
         LIEGE_L = WHERE(ABS(WAVELENGTH(*) - X_L) EQ MIN(ABS(WAVELENGTH(*) - X_L) ) )
         LIEGE_R = WHERE(ABS(WAVELENGTH(*) - X_R) EQ MIN(ABS(WAVELENGTH(*) - X_R) ) )

         ADU_INTENSITY = MEDIAN(REF_spectra(DATA_L(0):DATA_R(0),LINE_REF))

         IF PROFILE_METHOD EQ 1 THEN BEGIN ; CALCULATION USING LIEGE PROFILE
          LIEGE_INTENSITY = MEDIAN(DARKENING_FACTOR*INTENSITY(LIEGE_L(0):LIEGE_R(0)))
          FACTOR_CALCUL(K) = LIEGE_INTENSITY / ADU_INTENSITY
         ENDIF

         IF PROFILE_METHOD EQ 2 THEN BEGIN ; CALCULATION USING DAVID PROFILE CORRECTED FOR LIMB DARKENING
          INTENSITY_REF_SPLINE = spline(WAVELENGTH_REF,INTENSITY_REF,WAVELENGTH_DATA(*,LINE_REF),10.)
          ; LIMITS OF AREAS:
          DAVID_L = WHERE(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_L) EQ MIN(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_L) ) )
          DAVID_R = WHERE(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_R) EQ MIN(ABS(WAVELENGTH_DATA(*,LINE_REF) - X_R) ) )

          DAVID_INTENSITY = MEDIAN(INTENSITY_REF_SPLINE(DAVID_L(0):DAVID_R(0)))
          FACTOR_CALCUL(K) = DAVID_INTENSITY / ADU_INTENSITY
         ENDIF

       ENDFOR


       CASE NO_AREA OF
         1:       BEGIN
                 FACTOR_DATA_TEMP = FACTOR_CALCUL
                 FAC_DEGREE = 0
                 FACTOR_DATA = FLTARR(3)
                 FACTOR_DATA(0) = FACTOR_DATA_TEMP(0)
                 FACTOR_DATA(1) = 0.
                 FACTOR_DATA(2) = 0.
              END
         2:       BEGIN
                 FACTOR_DATA_TEMP = POLY_FIT(DATA_C,FACTOR_CALCUL,1,/DOUBLE,YFIT=FACTOR_POLY)
                 FAC_DEGREE = 1
                 FACTOR_DATA = FLTARR(3)
                 FACTOR_DATA(0) = FACTOR_DATA_TEMP(0)
                 FACTOR_DATA(1) = FACTOR_DATA_TEMP(1)
                 FACTOR_DATA(2) = 0.
              END
         3:       BEGIN
                 FACTOR_DATA = POLY_FIT(DATA_C,FACTOR_CALCUL,2,/DOUBLE,YFIT=FACTOR_POLY)
                 FAC_DEGREE = 2
              END
         ELSE:    BEGIN
                 FACTOR_DATA = POLY_FIT(DATA_C,FACTOR_CALCUL,2,/DOUBLE,YFIT=FACTOR_POLY)
                 FAC_DEGREE = 2
              END
       ENDCASE

       X_COORD = FINDGEN(1280)
       FACTOR_DATA_VECTOR(*,LINE_REF) = 0.
       FACTOR_DATA_COEF(*,LINE_REF) = 0.

       FOR H=0,FAC_DEGREE DO BEGIN
         FACTOR_DATA_VECTOR(*,LINE_REF) = FACTOR_DATA_VECTOR(*,LINE_REF) + X_COORD^H*FACTOR_DATA(H)
         FACTOR_DATA_COEF(H,LINE_REF) = FACTOR_DATA(H)
       ENDFOR

       ;INTENSITY CORRECTION OF SPECTRUM:
       ;FACTOR_DATA(LINE_REF) = 1./FACTOR ;TAK MA BYC
       Y_TITLE(LINE_REF) = 'PERCENT OF SOLAR CONTINUUM'

       WAIT, 0.5

       REF_EXP_TIME(LINE_REF) = EXP_TIME(LINE_REF)
       REF_TIME_CONTROL(LINE_REF) = 1

         HEADER_TEMP = HEADER_SPECTRA(*,LINE_REF)
                                                    ;45678901234567890123456789012345678901234567890
         FXADDPAR,HEADER_TEMP,'LAMBDA 0',DISPERSION_COEF(0,LINE_REF),'wavelength of the first pixel [Angstroems] MZ'
         FXADDPAR,HEADER_TEMP,'DISPCOEF',DISPERSION_COEF(1,LINE_REF),'dispert curve linear coeff [Angstroems/pxl] MZ'
         ;FXADDPAR,HEADER_TEMP,'DEG INT',FAC_DEGREE,              'degree of poly fit of intensity calib coeff MZ'
         FOR G=0,2 DO BEGIN;45678901234567890123456789012345678901234567890
          IF G EQ 0 THEN COMMENT='calib coef [% sol cont/ADU] of the 1st pxl MZ'
          IF G EQ 1 THEN COMMENT='calib coef - linear change MZ'
          IF G EQ 2 THEN COMMENT='calib coef - quadratic change MZ'
          FXADDPAR,HEADER_TEMP,'INT'+STRCOMPRESS(STRING(G),/REMOVE_ALL)+'COEF',FACTOR_DATA(G)*FACTOR_EXP(LINE_REF),COMMENT
         ENDFOR
         writefits,data_name(LINE_REF),data(*,*,LINE_REF),HEADER_TEMP ; modification HEADERS OF SPECTRAL IMAGES


       ;MENU,['YES','NO'],QQQ,TITLE='MAKE THIS IMAGE REFERENCE IMAGE?',BUTSIZE=300,xoffset=0,yoffset=-15
       ;IF QQQ EQ 1 THEN BEGIN
       ; OPENW,1,'calibration_coefficients.txt'
       ; FOR I=0,NO_CAM-1 DO FOR J=0,1 DO PRINTF,1,DISPERSION_COEF(J,I)
       ; FOR I=0,NO_CAM-1 DO PRINTF,1,FACTOR_DATA(I)
       ; FOR I=0,NO_CAM-1 DO PRINTF,1,REF_EXP_TIME(I)
       ; CLOSE,1
       ;ENDIF


       CALIBRATION_COMPLETE = 1



       LOADCT, 3





    ENDIF

    IF (LETR EQ 8) THEN BEGIN
       MENU,[STRUPCASE(CAM_INDEX[1:NO_CAM-1])],PROF_CAM,TITLE='CHOOSE CAMERA?',BUTSIZE=240,xoffset=0,yoffset=-15
       WSET,PROF_CAM & WSHOW,PROF_CAM
       PROFILES, trans_image_cut(*,*,PROF_CAM)
    ENDIF

    IF (LETR EQ 9) THEN BEGIN

    IF FIRST_AGAIN EQ 0 THEN BEGIN
       MENU,['OK'],PROF_CAM,TITLE='YOU HAVE TO EXTRACT SPECTRUM TO CALIBRATION',BUTSIZE=600,xoffset=0,yoffset=-15
       GOTO,START
    ENDIF

       CLOSE,/ALL
       ON_IOERROR, A223
       A1=' '
       OPENR,1,'liege.txt'
       FOR I=0.,999999. DO READF,1,A1
       A223: no_lines=I
       CLOSE,1
       ON_IOERROR,NULL

       WHOLE_SPECTRUM = DBLARR(2,NO_LINES)
       OPENR,1,'liege.txt'
       READF,1,WHOLE_SPECTRUM
       CLOSE,1

       WAVELENGTH = WHOLE_SPECTRUM(0,*)
       INTENSITY  = WHOLE_SPECTRUM(1,*)

       DISPERSION_COMPLETE = 0

       ;MENU,['HA','CA','D3','HB','IR'],LINE_CORR,TITLE='CHOOSE CAMERA?',BUTSIZE=300,xoffset=0,yoffset=-15
       MENU,[STRUPCASE(CAM_INDEX[1:NO_CAM-1])],LINE_CORR,TITLE='CHOOSE CAMERA?',BUTSIZE=300,xoffset=0,yoffset=-15

       CASE CAM_INDEX(LINE_CORR) OF
         'ha': LINE_CORR_LIEGE = 1
         'ca': LINE_CORR_LIEGE = 2
         'd3': LINE_CORR_LIEGE = 3
         'hb': LINE_CORR_LIEGE = 4
         'ir': LINE_CORR_LIEGE = 5
       ENDCASE

		IF CAM_INDEX(LINE_CORR) EQ 'ca' THEN BEGIN

			menu,['K-3933.682','H-3968.492'],LIN,TITLE='CHOOSE CALCIUM LINE',BUTSIZE=600
			IF LIN EQ 1 THEN BEGIN
				;for Ca K line
				MIN_WL = [0.,6553.,3919.,5863.,4848.,8529.5]
				MAX_WL = [0.,6573.,3946.,5888.,4875.,8555.5]
				LAMBDA_CENTR = [0.0, 6562.808, 3933.682, 5875.618, 4861.342, 8542.144  ]
			ENDIF

			IF LIN EQ 2 THEN BEGIN
				;for Ca H line
				MIN_WL = [0.,6553.,3955.,5863.,4848.,8529.5]
				MAX_WL = [0.,6573.,3982.,5888.,4875.,8555.5]
				LAMBDA_CENTR = [0.0, 6562.808, 3968.492, 5875.618, 4861.342, 8542.144  ]
			ENDIF

		ENDIF


       REPEAT BEGIN

         INDEX_WL_MIN=WHERE(WAVELENGTH EQ MIN_WL(LINE_CORR_LIEGE))  ;LIEGE SPECTRUM
         INDEX_WL_MAX=WHERE(WAVELENGTH EQ MAX_WL(LINE_CORR_LIEGE))

         window,NO_CAM+1,ypos=502,xpos=80,xsize=1280,ysize=500,title=' REFERENCE LIEGE SPECTRUM --  CLICK ON CHOSEN LINE '
; PLOTTING IN RANGE THAT ANNOTATIONS OF REFERENCE LINES WILL BE VISIBLE:
         IF LINE_CORR EQ 2 THEN BEGIN ; CA LINE
          PLOT,WAVELENGTH[INDEX_WL_MIN(0):INDEX_WL_MAX(0)],INTENSITY[INDEX_WL_MIN(0):INDEX_WL_MAX(0)],$
                   XRANGE=[WAVELENGTH(INDEX_WL_MIN(0)),WAVELENGTH(INDEX_WL_MAX(0))],XSTYLE=1, YRANGE=[-50,100]
         ENDIF ELSE BEGIN
          PLOT,WAVELENGTH[INDEX_WL_MIN(0):INDEX_WL_MAX(0)],INTENSITY[INDEX_WL_MIN(0):INDEX_WL_MAX(0)],$
                   XRANGE=[WAVELENGTH(INDEX_WL_MIN(0)),WAVELENGTH(INDEX_WL_MAX(0))],XSTYLE=1
         ENDELSE
         for j=0,NO_REF_LINES-1 do begin
          XYOUTS,lines(J).wave,lines(j).inten-2.,lines(j).name+'___',/DATA,ORIENTATION=90,ALIGNMENT=1.0,color=lines(j).atm_color
         endfor

         IF DISPERSION_COMPLETE EQ 0 THEN BEGIN
          window,no_cam,ypos=0,xpos=80,xsize=1280,ysize=500,title=' DATA SPECTRUM -- CLICK ON CHOSEN LINE '
          PLOT,spectra(*,LINE_CORR), XRANGE=[0,1280], XSTYLE=1, TITLE=strupcase(CAM_INDEX(LINE_CORR))+' SPECTRUM ',charsize=1.0
         ENDIF

         IF DISPERSION_COMPLETE EQ 1 THEN BEGIN
          WINDOW,12,ypos=502,xpos=80,xsize=500,ysize=500,TITLE='DISPERSION CURVE'
          PLOT,X_DATA,LINES(LINE_INDEX).WAVE,/YNOZERO,PSYM=1,SYMSIZE=1.5,XRANGE=[0,1280],XSTYLE=1
          OPLOT,INDGEN(1280)*DISPERSION_COEF(1,LINE_CORR)+DISPERSION_COEF(0,LINE_CORR),COLOR=150
          XYOUTS,0.18,0.90, 'LAMBDA(0)  =' +STRCOMPRESS(STRING(DISPERSION_COEF(0,LINE_CORR),format='(F10.4)'))+' +/-'+STRCOMPRESS(STRING(ERROR(0),format='(F10.4)'))+' A',/NORMAL,COLOR=150
          XYOUTS,0.18,0.85,'DISPERSION ='+STRCOMPRESS(STRING(DISPERSION_COEF(1,LINE_CORR),format='(F10.7)'))+' +/-'+STRCOMPRESS(STRING(ERROR(1),format='(F10.7)'))+' A/PXL',/NORMAL,COLOR=150
          XYOUTS,0.18,0.80,'WHOLE WIDTH ERROR =' + STRCOMPRESS(STRING(ERROR(1)*1280.,format='(F10.7)'))+' A ',/NORMAL,COLOR=150
          XYOUTS,0.18,0.75,'WHOLE WIDTH ERROR =' + STRCOMPRESS(STRING((ERROR(1)*1280.+ERROR(0))/DISPERSION_COEF(1,LINE_CORR),format='(F10.7)'))+' PXL ',/NORMAL,COLOR=150
          PRINT,'LAMBDA(0) =' +STRCOMPRESS(STRING(DISPERSION_COEF(0,LINE_CORR),format='(F10.4)'))+' +/-'+STRCOMPRESS(STRING(ERROR(0),format='(F10.4)'))+' A'
          PRINT,'DISPERSION ='+STRCOMPRESS(STRING(DISPERSION_COEF(1,LINE_CORR),format='(F10.7)'))+' +/-'+STRCOMPRESS(STRING(ERROR(1),format='(F10.7)'))+' A/PXL'
          window,no_cam,ypos=0,xpos=80,xsize=1280,ysize=500,title=' DATA SPECTRUM -- AFTER DISPERSION '
          PLOT,WAVELENGTH_DATA(*,LINE_CORR),spectra(*,LINE_CORR),XRANGE=[DISPERSION_COEF(0,LINE_CORR),DISPERSION_COEF(1,LINE_CORR)*1279.+DISPERSION_COEF(0,LINE_CORR)], $
              XSTYLE=1, TITLE=strupcase(CAM_INDEX(LINE_CORR))+' SPECTRUM ',charsize=1.0
         ENDIF

;nie wyswietla sie dziadostwo

         MENU,['DISPERSION CURVE','REVERSE','OK'],CONTROL,TITLE='AGAIN?',BUTSIZE=160,xoffset=0,yoffset=0

         IF CONTROL EQ 1 THEN BEGIN
          MENU,['1','2','3','4','5'],NO_LI,TITLE='HOW MANY LINES?',BUTSIZE=160,xoffset=0,yoffset=0
          LINE_INDEX = INTARR(NO_LI)
          X_DATA = FLTARR(NO_LI)
          X_DATA_POLY = FLTARR(NO_LI)
          Y_DATA = FLTARR(NO_LI)
          FOR I=0,NO_LI-1 DO BEGIN
              window,NO_CAM+1,ypos=502,xpos=80,xsize=1280,ysize=500,title=' REFERENCE LIEGE SPECTRUM --  CLICK ON CHOSEN LINE '
              PLOT,WAVELENGTH[INDEX_WL_MIN(0):INDEX_WL_MAX(0)],INTENSITY[INDEX_WL_MIN(0):INDEX_WL_MAX(0)],$
                             XRANGE=[WAVELENGTH(INDEX_WL_MIN(0)),WAVELENGTH(INDEX_WL_MAX(0))],XSTYLE=1,YRANGE=[-50,100]

              for j=0,NO_REF_LINES-1 do begin ; LINE INFORMATION
                 XYOUTS,lines(J).wave,lines(j).inten-2.,lines(j).name+'___',/DATA,ORIENTATION=90,ALIGNMENT=1.0,color=lines(j).atm_color
              endfor

              FOR R=0,I-1 DO BEGIN ;CHOSEM LINES INFORMATION
                 XYOUTS,lines(LINE_INDEX(R)).wave,lines(LINE_INDEX(R)).inten-2.,'|___________________________|',$
                    /DATA,ORIENTATION=90,ALIGNMENT=1.0,color=lines(LINE_INDEX(R)).atm_color
              ENDFOR

              CURSOR,X_REF,Y_REF,/DATA,/DOWN ;CLICK ON REFERENCE LINE ONLY WITH LIMITED ACCURACY
              PRINT,X_REF
              DIFF=FLTARR(NO_REF_LINES) ;TOTAL NO LINES IN CHOSEN CAMERA
              FOR J=0,NO_REF_LINES-1 DO BEGIN
                 DIFF(J) = ABS(X_REF - lines(j).wave)
              ENDFOR
              PRINT,DIFF
              LINE_INDEX(I) = WHERE(DIFF EQ MIN(DIFF))
              PRINT,LINE_INDEX
              XYOUTS,lines(LINE_INDEX(I)).wave,lines(LINE_INDEX(I)).inten-2.,'|___________________________|',$
                    /DATA,ORIENTATION=90,ALIGNMENT=1.0,color=lines(LINE_INDEX(I)).atm_color


              window,no_cam,ypos=0,xpos=80,xsize=1280,ysize=500,title=' DATA SPECTRUM -- CLICK ON CHOSEN LINE '
              PLOT,spectra(*,LINE_CORR), XRANGE=[0,1280], XSTYLE=1, TITLE=strupcase(CAM_INDEX(LINE_CORR))+' SPECTRUM ',charsize=1.0

              FOR R=0,I-1 DO BEGIN
                 PLOTS,[x_DATA(R),x_DATA(R)],[0,Y_DATA(R)],/DATA,color=lines(LINE_INDEX(R)).atm_color
              ENDFOR

              CURSOR,X_,Y_,/DATA,/DOWN
              X_DATA(I)=X_
              Y_DATA(I)=Y_
              PLOTS,[x_DATA(I),x_DATA(I)],[0,Y_DATA(I)],/DATA,color=lines(LINE_INDEX(I)).atm_color

              ;POLY FITTING TO CHOSEN LINES:
              Y_TEMP = INDGEN (1280)
              RANGE = 6
              X_MIN_POS = WHERE( SPECTRA[ X_DATA(I)-RANGE:X_DATA(I)+RANGE, LINE_CORR ] EQ MIN ( SPECTRA[ X_DATA(I)-RANGE:X_DATA(I)+RANGE, LINE_CORR ] ) )
              X_MIN_POS = X_MIN_POS(0) + X_DATA(I)-RANGE
              TEMP = POLY_FIT(Y_TEMP[X_MIN_POS-2:X_MIN_POS+2],spectra(X_MIN_POS-2:X_MIN_POS+2,LINE_CORR),2, /DOUBLE)
              X_DATA_POLY(I) =  -TEMP(1)/(2*TEMP(2))
              ;STOP


          ENDFOR

          WAIT,0.5
         ;STOP
          DISPERSION_COEF(*,LINE_CORR) = POLY_FIT(X_DATA_POLY,LINES(LINE_INDEX).WAVE,1,/DOUBLE,SIGMA=ERROR) ;LINEAR DISPERSION ???

          DISPERSION_COMPLETE = 1

          WAVELENGTH_DATA(*,LINE_CORR)=FINDGEN(1280.)*DISPERSION_COEF(1,LINE_CORR) + DISPERSION_COEF(0,LINE_CORR)

         ENDIF


         IF CONTROL EQ 2 THEN BEGIN

          spectra(*,LINE_CORR)=REVERSE(SPECTRA(*,LINE_CORR)); MAYBE BETTER INDIVIDUALLY FOR EACH LINE ??
          REVERSE_CONTROL(LINE_CORR) = REVERSE_CONTROL(LINE_CORR) + 1
          ;window,no_cam,ypos=80,xpos=80,xsize=320*4,ysize=128*4,title=' DATA SPECTRUM -- CLICK ON CHOSEN LINE '
          ;PLOT,spectra(*,LINE_CORR), XRANGE=[0,1280], XSTYLE=1, TITLE=strupcase(CAM_INDEX(LINE_CORR))+' SPECTRUM ',charsize=2.0
          ;IF REV(LINE_CORR) EQ 0 THEN REV(LINE_CORR) = 1
          ;IF REV(LINE_CORR) EQ 1 THEN REV(LINE_CORR) = 0

         ENDIF

       ENDREP UNTIL CONTROL EQ 3
    ENDIF


    IF (LETR EQ 10) THEN BEGIN ;SCATTTERED LIGHT SUBTRACTION

              MENU,['SCAT LIGHT MATRIX','CONTINUUM LEVEL'],$
              SCAT_METHOD,TITLE='CHOOSE METHOD',BUTSIZE=300,xoffset=0,yoffset=0


                 REF_SCATT_LIGHT = FLTARR(1280, NO_CAM)

                 ;calculation no_pts
                 CLOSE,/ALL
                 ON_IOERROR, D888
                 AAA=' '
                 OPENR,5,'scat_light_position_'+strlowcase(CAM_INDEX(1))+'.txt'
                 FOR I=0.,999. DO READF,5,AAA
                 D888: NO_PTS=I
                 CLOSE,5
                 ON_IOERROR,NULL

                 CORR_SCATT_LIGHT = DBLARR(1280, NO_CAM, NO_PTS)
                 SCAT_LIGHT_SPECTRA = DBLARR(1280, NO_CAM, NO_PTS)
                 SCAT_LIGHT_POSITION = DBLARR(3,NO_CAM,NO_PTS)
                 CONTINUUM_LIMITS = INTARR(2,NO_CAM)

                 SCAT_NAME = FINDFILE('scat_light_position.txt')


                 ;IF SCAT_NAME(0) NE '' THEN BEGIN

                   FOR NN=1,NO_CAM-1 DO BEGIN

                    WHOLE_FILE = DBLARR(1280,NO_PTS)

                    OPENR,4,'scat_light_profiles_'+strlowcase(CAM_INDEX(NN))+'.txt'
                    READF,4,WHOLE_FILE
                    CLOSE,4

                    SCAT_LIGHT_SPECTRA(*, NN, *) =  WHOLE_FILE(*,*); (ROW, CAM, POINT)

                    WHOLE_FILE = DBLARR(3,NO_PTS)

                    OPENR,5,'scat_light_position_'+strlowcase(CAM_INDEX(NN))+'.txt'
                    READF,5,WHOLE_FILE  ;; NOW HERE ARE STORED POSITIONS OF SAMPLED SCATTERED LIGHT
                    CLOSE,5

                   ENDFOR

                 ;ENDIF ;ELSE BEGIN
                   ;MENU,['OK'],PROF_CAM,TITLE='YOU HAVE TO CALCULATE SCATTERED LIGHT TO SUBTRACT IT',BUTSIZE=800,xoffset=0,yoffset=-15
                   ;GOTO,AFTER_SCAT_LIGHT_REMOVED
                 ;ENDELSE
              ;STOP

;uuuuuuuuuuuuuuuuuuuuuuuuuuu          tutaj          uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu

              IF SCAT_METHOD EQ 1 THEN BEGIN

                     SCAT_PARAM = replicate({fixed:0, limited:[0,0], limits:[0.D,0.D]},2)
                 SCAT_MATRIX = FLTARR(1280,512,NO_CAM)
                 SCAT_COEFF_MATRIX = FLTARR(2,1280)

                 POS_SCAT  = WHOLE_FILE(1,*)
                 DISTANCE_R_SCAT = SQRT( (SUN_X_CENTR-X_SLIT_POSITION)^2. + (SUN_Y_CENTR-POS_SCAT*2.)^2. ) / SUN_RADIUS

                 xx=(findgen(100)+1.)/1000.+1.
                 X_DIST = FLTARR(512)
                 FOR I = 0, 511 DO BEGIN
                   X_DIST(I) = SQRT( (SUN_X_CENTR-X_SLIT_POSITION)^2. + (SUN_Y_CENTR-I*2.)^2. ) / SUN_RADIUS
                 ENDFOR

                 FOR JJ=1,NO_CAM-1 DO BEGIN

                   ;CONTOL IF FILE WITH SCAT MATRIX EXIST
                   FILE_CONTROL = FINDFILE('resources\scat_light_matrix_'+data_time+'_'+cam_index(jj)+'.fit')

                   ;READING FILE WITH SCAT MATRIX
                   IF FILE_CONTROL(0) EQ 'resources\scat_light_matrix_'+data_time+'_'+cam_index(jj)+'.fit' THEN BEGIN
                    SCAT_MATRIX(*,*,JJ)=readfits('resources\scat_light_matrix_'+data_time+'_'+cam_index(jj)+'.fit',HH)
                   ENDIF ELSE BEGIN

                    window,10, TITLE='CALCULATING SCATTERED LIGHT MATRIX'
                    FOR I=0,1279 DO BEGIN

                        VERT_SCAT = SCAT_LIGHT_SPECTRA (I,JJ,*) ;/ MEDIAN(SCAT_LIGHT_SPECTRA (I,JJ,*))
                        COL_RANGE = MAX(SCAT_LIGHT_SPECTRA(*,JJ,*))

                        ;CONTROLLING LIMITS OF PARAMETERS
                        ;SEE: http://www.physics.wisc.edu/~craigm/idl/mpfittut.html

                        SCAT_PARAM(0).LIMITED(0) = 0
                        ;SCAT_PARAM(0).LIMITS(0) =  -2.1

                        SCAT_PARAM(0).LIMITED(0) = 0
                        ;SCAT_PARAM(0).LIMITS(1) =  -2.4

                        SCAT_PARAM(1).LIMITED(0) = 0
                        ;SCAT_PARAM(1).LIMITS(0) = ALOG (2.6)

                        SCAT_PARAM(1).LIMITED(0) = 0
                        ;SCAT_PARAM(1).LIMITS(1) = ALOG (3.4)
                        scat_param(1).fixed = 0

                        START = [150.,-0.4]

                        SCAT_COEFF = MPFITFUN('SCAT_LOG', DISTANCE_R_SCAT[0,*], VERT_SCAT[0,0,*] , 1.0, START, PARINFO=SCAT_PARAM)
                        SCAT_COEFF_MATRIX(*,I) = SCAT_COEFF

                        ;X_DIST = X_DIST ; ARBITRARY CUTTING LEVEL TO PREVENT FLOATING OVERFLOW DURING ALOG
                        COL = SCAT_COEFF(0) * (SCAT_COEFF(1) - ALOG10(X_DIST>1.001-1))
                        SCAT_MATRIX(I,*,JJ) = COL
                        PLOT, X_DIST, COL,YRANGE=[0,1.5*COL_RANGE],XTITLE='DISTANCE [SOLAR RADII]', YTITLE='SCATTERED LIGHT LEVEL', $
                             TITLE=CAM_INDEX(JJ)+'  LAMBDA'+STRCOMPRESS(STRING(WAVELENGTH_DATA[I,JJ]))+STRING(I)+' / 1280',XSTYLE=1
                        OPLOT,DISTANCE_R_SCAT,VERT_SCAT,PSYM=2,COLOR=150
                        ;WAIT,0.0

                    ENDFOR
                    ;STOP
                    ;

                   ENDELSE

                   window,0,ypos=0*150,xpos=0*75,xsize=320*4,ysize=128*4,title=strUPcase(CAM_INDEX(JJ))+' SCATTERED LIGHT MATRIX'
                   TVSCL , SCAT_MATRIX(*,*,JJ)
                   MENU,['OK','NO','CANCEL'],$
                   SSS,TITLE='USE THIS MATRIX TO SUBTRACTION?',BUTSIZE=400,xoffset=0,yoffset=0
                   IF SSS EQ 1 THEN BEGIN
                    TRANS_IMAGE_CUT(*,*,JJ)    = TRANS_IMAGE_CUT(*,*,JJ)      - SCAT_MATRIX(*,*,JJ)
                    writefits, 'resources\scat_light_matrix_'+data_time+'_'+STRlowCASE(cam_index(jj))+'.fit', SCAT_MATRIX(*,*,JJ)
                   ENDIF
                   IF SSS EQ 3 THEN BREAK
                 ENDFOR

              GOTO,AFTER_SCAT_LIGHT_REMOVED
              ENDIF


;uuuuuuuuuuuuuuuuuuuuuuuuuuu          tutaj          uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu


          IF SCAT_METHOD EQ 2 THEN BEGIN ;COMPARE METHOD

                 IF SCAT_NAME(0) EQ 'scat_light_position.txt' THEN BEGIN

                    OPENR, 2, 'scat_light_position.txt'
                    READF, 2, CONTINUUM_LIMITS
                    CLOSE,2

                 ENDIF ELSE BEGIN

                    FOR NN=1,NO_CAM-1 DO BEGIN
                        SCAT_LIGHT_POSITION(*,NN,*) = WHOLE_FILE(*,*) ; ([X,Y,DISTANCE],CAM_INDEX,POINTS_INDEX)
                        window,no_cam+NN,ypos=80,xpos=80,xsize=320*4,ysize=128*7,title=' SPECTRUM OF SCATTERED LIGHT FROM MEMORY - MARK CONTINUUM LIMITS'
                        PLOT, SMOOTH(SCAT_LIGHT_SPECTRA(*,NN,0),11), PSYM = 3, XSTYLE=1, YRANGE = [0,MAX(SCAT_LIGHT_SPECTRA(*,NN,*))]
                        FOR J=1, NO_PTS-1 DO BEGIN
                           OPLOT, SMOOTH(SCAT_LIGHT_SPECTRA(*,NN,J),11), PSYM =3, COLOR=J*(256./NO_PTS-1)
                        ENDFOR

                        CURSOR, X_1,Y,/DATA,/DOWN
                        PLOTS,[X_1,X_1],[0,Y],/DATA,COLOR=150

                        CURSOR, X_2,Y,/DATA,/DOWN
                        PLOTS,[X_2,X_2],[0,Y],/DATA,COLOR=150
                        WAIT,0.4

                        X_LEFT =  X_1 < X_2
                        X_RIGHT = X_1 > X_2

                        X_CONT_LEFT = ROUND(X_LEFT)
                        X_CONT_RIGHT = ROUND(X_RIGHT)
                        CONTINUUM_LIMITS(0,NN)=X_CONT_LEFT
                        CONTINUUM_LIMITS(1,NN)=X_CONT_RIGHT
                    ENDFOR

                    OPENW, 2, 'scat_light_position.txt'
                    PRINTF, 2, CONTINUUM_LIMITS
                    CLOSE, 2

                 ENDELSE



                    ;CALCULATION OF CORRECTION FACTOR OF CONTINUUM LEVEL
                    CONT_FACTOR = DBLARR(NO_PTS)

                 FOR NN=1,NO_CAM-1 DO BEGIN

                    FOR J=1, NO_PTS-1 DO BEGIN
                        CONT_FACTOR(J) = MEDIAN(SCAT_LIGHT_SPECTRA(CONTINUUM_LIMITS(0,NN):CONTINUUM_LIMITS(1,NN),NN,J))
                        CORR_SCATT_LIGHT(*,NN,J) = SCAT_LIGHT_SPECTRA(*,NN,J) / CONT_FACTOR(J)
                    ENDFOR

                    REF_SCATT_LIGHT(*, NN) = TOTAL(CORR_SCATT_LIGHT(*,NN,*),3 ) / NO_PTS

                    FOR ROW=0,511 DO BEGIN
                        CURRENT_FACTOR = MEDIAN(TRANS_IMAGE_CUT(CONTINUUM_LIMITS(0,NN):CONTINUUM_LIMITS(1,NN),ROW,NN))
                        TRANS_IMAGE_CUT(*,ROW,NN)     = TRANS_IMAGE_CUT(*,ROW,NN)    - REF_SCATT_LIGHT(*,NN)*CURRENT_FACTOR
                        TEMP_FOR_INTEGRATION(*,ROW,NN) = TEMP_FOR_INTEGRATION(*,ROW,NN) - REF_SCATT_LIGHT(*,NN)*CURRENT_FACTOR
                    ENDFOR

                    window,NN,ypos=0*150,xpos=NN*75,xsize=320*4,ysize=128*4,title=strUPcase(CAM_INDEX(NN))+' SCATTERED LIGHT REMOVED ( '+data_name(0)+' )'
                    TVSCL , TRANS_IMAGE_CUT(*,*,NN)

                 ENDFOR

              FIRST_AGAIN = 1
              GOTO,AFTER_SCAT_LIGHT_REMOVED

          ENDIF
    ENDIF

    IF (LETR EQ 11) THEN BEGIN

       IF FIRST_AGAIN EQ 0 THEN BEGIN
         MENU,['OK'],PROF_CAM,TITLE='YOU HAVE TO EXTRACT SPECTRUM TO DOPPLER MEASUREMENTS',BUTSIZE=800,xoffset=0,yoffset=-15
         GOTO,START
       ENDIF

       VELOCITY = DBLARR(NO_CAM)

       REPEAT BEGIN

         MENU,['BISECTRIES','GAUSS FIT','OK AND SAVE','CANCEL'],DOPP_METHOD,TITLE='DOPPLER MEASUREMENTS',BUTSIZE=400,xoffset=0,yoffset=-15

         !P.MULTI = [0,2,3,0,0]
         window,no_cam,ypos=110,xpos=80,xsize=320*4,ysize=128*7,title=' DOPPLER VELOCITY MEASUREMENTS'

         LOADCT,39

         IF DOPP_METHOD EQ 2 THEN BEGIN  ;GAUSS FIT

          FOR LINE_DOPP=1,NO_CAM-1 DO BEGIN

              CASE CAM_INDEX(LINE_DOPP) OF
                 'ha': LINE_DOPP_MOD = 1
                 'ca': LINE_DOPP_MOD = 2
                 'd3': LINE_DOPP_MOD = 3
                 'hb': LINE_DOPP_MOD = 4
                 'ir': LINE_DOPP_MOD = 5
              ENDCASE

              R = GAUSSFIT( WAVELENGTH_DATA(*,LINE_DOPP), $
                 FACTOR_EXP(LINE_DOPP)*FACTOR_CGS(LINE_DOPP)*FACTOR_DATA_VECTOR(*,LINE_DOPP)*spectra(*,LINE_DOPP), $
                 GAUSS_COEF, NTERMS=4)

              VELOCITY(LINE_DOPP) =  -1. * (LAMBDA_CENTR(LINE_DOPP_MOD) - GAUSS_COEF(1)) /  LAMBDA_CENTR (LINE_DOPP_MOD) * 299792.458

              PLOT, WAVELENGTH_DATA(*,LINE_DOPP), $
                 FACTOR_EXP(LINE_DOPP)*FACTOR_CGS(LINE_DOPP)*FACTOR_DATA_VECTOR(*,LINE_DOPP)*spectra(*,LINE_DOPP),$
                 TITLE=STRUPCASE(CAM_INDEX(LINE_DOPP))+STRCOMPRESS(STRING(VELOCITY(LINE_DOPP)))+' KM/S',CHARSIZE=1.8, XSTYLE=1
              OPLOT, WAVELENGTH_DATA(*,LINE_DOPP), R, COLOR = 150

          ENDFOR

         ENDIF

         IF DOPP_METHOD EQ 1 THEN BEGIN ; BISECTRIES

          FOR LINE_DOPP=1,NO_CAM-1 DO BEGIN

              CASE CAM_INDEX(LINE_DOPP) OF
                 'ha': LINE_DOPP_MOD = 1
                 'ca': LINE_DOPP_MOD = 2
                 'd3': LINE_DOPP_MOD = 3
                 'hb': LINE_DOPP_MOD = 4
                 'ir': LINE_DOPP_MOD = 5
              ENDCASE

              V=FLTARR(4)
              VELOC_MODIFY,WAVELENGTH_DATA(*,LINE_DOPP),SPECTRA(*,LINE_DOPP),V,VELOCITY_TEMP,1,NO_CAM,LINE_DOPP_MOD
              VELOCITY(LINE_DOPP) = VELOCITY_TEMP
              ;PRINT,V_DOPP

          ENDFOR


         ENDIF

       ENDREP UNTIL DOPP_METHOD EQ 3 OR DOPP_METHOD EQ 4

       IF DOPP_METHOD EQ 3 THEN BEGIN
         ;format of output data,line by line:

         ;time (HHMMSS)
         ;time (seconds from the begin of the day)
         ;slit position (pxl)
         ;centrum of analyzed area (pxl)
         ;H-alpha calculated velocity
         ;Ca  calculated velocity
         ;D3  calculated velocity
         ;H-beta calculated velocity
         ;Ca IR   calculated velocity

         OPENW,1,'velocity.dat',/APPEND, WIDTH=160
         TIME=STRMID(data_time(0),0,2)+STRMID(data_time(0),3,2)+STRMID(data_time(0),6,2)

         HHORA=(STRMID(TIME,0,2))
         MMIN=(STRMID(TIME,2,2))
         SSEK=(STRMID(TIME,4,2))

         TOTSEK=HHORA*60.0*60.0+MMIN*60.0+SSEK

         PRINTF, 1, TIME, TOTSEK, X_SLIT_POSITION, (Y_UP+Y_DOWN)/2., VELOCITY[1:NO_CAM-1]
         CLOSE,1
       ENDIF


       !P.MULTI = 0
       LOADCT,3

    ENDIF


    IF (LETR EQ 12) THEN BEGIN ;SAVE


       IF SUN_X_CENTR NE -1. AND SUN_Y_CENTR NE -1. THEN BEGIN ; CALCULATE ONLY IF DISC CENTER IS CALCULATED

         ;MENU,['HA','CA','D3','HB','IR'],LINE_CORR,TITLE='CHOOSE CAMERA?',BUTSIZE=300,xoffset=0,yoffset=-15
         OPENW,1,'calibration_coefficients.txt'
         ;PRINTF,1,SLIT_WIDTH
         FOR I=0,NO_CAM-1 DO FOR J=0,1 DO PRINTF,1,DISPERSION_COEF(J,I)
         FOR I=0,NO_CAM-1 DO FOR H=0,2 DO PRINTF,1,FACTOR_DATA_COEF(H,I)
         FOR I=0,NO_CAM-1 DO PRINTF,1,REF_EXP_TIME(I)
         CLOSE,1

         ;writing into the header of SJ image, coefficients of position of calculated disc center
         FXADDPAR,HEADER_SJ,'DISCX',SUN_X_CENTR,'x sun center in ref.frame of the image [pxl] MZ'
         FXADDPAR,HEADER_SJ,'DISCY',SUN_Y_CENTR,'y sun center in ref.frame of the image [pxl] MZ'
         FXADDPAR,HEADER_SJ,'SLIT POS',X_SLIT_POSITION,'x slit position on the image [pxl] MZ'
         writefits,data_name(0),data(*,*,0),HEADER_SJ ; modyfication only in header in SJ image

       ENDIF ELSE BEGIN
         MENU,['OK'],LINE_CORR,TITLE='PLEASE CALCULATE DISC CENTER BEFORE SAVING COEFFICIENTS',BUTSIZE=600,xoffset=300,yoffset=300
       ENDELSE


    ENDIF

    IF (LETR EQ 13) THEN BEGIN  ;credits

         window,30,xpos=750,ypos=0,xsize=500,ysize=500,title='CREDITS'

         A=1.6
         B=1.5

         XYOUTS, 150*A, 290*A, 'HSFA-2 SOFTWARE PACKAGE ',                 CHARTHICK=1.5,  CHARSIZE=2.,    COLOR=255, ALIGNMENT=0.5,/device

         XYOUTS, 150*A, 245*A, 'author: ',                           CHARTHICK=1,    CHARSIZE=1.7,     COLOR=255, ALIGNMENT=0.5,/device

         XYOUTS, 150*A, 225*A, 'Maciej Zapior (1,2)',                         CHARTHICK=1,   CHARSIZE=1.7,    COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A, 205*A, '(1) - Astronomocal Institute of Wroclaw University',      CHARTHICK=1,   CHARSIZE=1.0,    COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A, 195*A, '(2) - Astronomocal Institute of the Academy of Science of the Czech Republic',$
                                                        CHARTHICK=1,     CHARSIZE=1.0,  COLOR=255, ALIGNMENT=0.5,/device

         XYOUTS, 150*A, 150*A, 'Thanx: ',                              CHARTHICK=1,     CHARSIZE=1.7,  COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A, 130*A, 'Pavel Kotrc (2)',                            CHARTHICK=1,  CHARSIZE=1.5,   COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A, 110*A, 'Petr Heizel (2)',                            CHARTHICK=1,  CHARSIZE=1.5,   COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A,  90*A, 'Petr Skoda (2)',                           CHARTHICK=1,   CHARSIZE=1.5,    COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A,  70*A, 'Jan Rybak',                          CHARTHICK=1,    CHARSIZE=1.5,     COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A,  50*A, 'Ales Kucera',                            CHARTHICK=1,  CHARSIZE=1.5,   COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A,  30*A, 'Pawel Rudawy (1)',                          CHARTHICK=1,     CHARSIZE=1.5,  COLOR=255, ALIGNMENT=0.5,/device
         XYOUTS, 150*A,  10*A, 'Arek Berlicki (1)',                        CHARTHICK=1,    CHARSIZE=1.5,     COLOR=255, ALIGNMENT=0.5,/device

       MENU,['OK','PLAY SNAKE'],SNAKE_CONTROL,TITLE='CONTINUE?',butsize=250
       IF SNAKE_CONTROL EQ 2 THEN SNAKE


    ENDIF

    IF (LETR EQ 14) THEN BEGIN  ;CGS AND INTEGRAL



       INT_SIGNAL_SLIT = FLTARR(512,NO_CAM)


       MENU,[STRUPCASE(CAM_INDEX[1:NO_CAM-1])],J,TITLE='CHOOSE CAMERA?',BUTSIZE=300,xoffset=0,yoffset=0
       Y_TITLE(J) = 'erg * sr^(-1) * cm^(-2) * s^(-1) * A^(-1) '  ;  <<< ----------    SPRAWDZIC

       IF CGS_DONE(J) EQ 0 THEN BEGIN
         ;FACTOR_DATA(J) = FACTOR_DATA(J) * 4.077 * 10.^(-5) ; <<<<<<<<<<<<<<   -------------    CHECK UNITS!!!!!
         FACTOR_CGS(J) = FACTOR_CGS(J) * CONTINUUM(LAMBDA_CENTR(J)) * 10.^(-2.) ; <<<<<<<<<<<<<<   -------------    CHECK UNITS!!!!!
         ;                                 CGS            'CAUSE WE HAVE 100 AS SOLAR CONTINUUM
         ; HERE IT IS PER 1A
         ; IN FACT THIS IS NOT IN CGS, BECAUSE IT IS NOT PER 1Hz

         CGS_DONE(J) = 1
       ENDIF
       ;PRINT, FACTOR_DATA
       PRINT, FACTOR_CGS

       LIMITS_CONTROL = 0

       REPEAT BEGIN
         ;window,NO_CAM,ypos=80,xpos=80,xsize=320*4,ysize=128*7,title=strupcase(CAM_INDEX(J))+' SPECTRUM'
         ;PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA(J)*spectra(*,j), XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE

         MENU,['MARK INTEGRAL LIMITS AND CALCULATE ','OK','SMOOTH -3-','SMOOTH -5-','SMOOTH -7-','SMOOTH -9-',$
                 'CALCULATE FOR WHOLE SLIT','CORRECTION OF OFFSET','SAVE INTEGRAL'],CONTROL,TITLE=CGS_TITLE,BUTSIZE=400,xoffset=0,yoffset=0

         IF CONTROL EQ 1 THEN BEGIN

         window,NO_CAM,ypos=80,xpos=80,xsize=320*4,ysize=128*7,title=' MARK INTEGRAL LIMITS'
         PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), $
         XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
         PLOTS,[MIN_WL(J),MAX_WL(J)],[0,0],/DATA,COLOR=150

         CURSOR, X_1,Y,/DATA,/DOWN
         PLOTS,[X_1,X_1],[0,Y],/DATA,COLOR=150

         CURSOR, X_2,Y,/DATA,/DOWN
         PLOTS,[X_2,X_2],[0,Y],/DATA,COLOR=150
         WAIT,0.4

         X_LEFT =  X_1 < X_2
         X_RIGHT = X_1 > X_2

         X_INDEX_LEFT =  WHERE( ABS(WAVELENGTH_DATA(*,J)-X_LEFT)  EQ  MIN(ABS(X_LEFT -WAVELENGTH_DATA(*,J))))
         X_INDEX_RIGHT = WHERE( ABS(WAVELENGTH_DATA(*,J)-X_RIGHT) EQ  MIN(ABS(X_RIGHT-WAVELENGTH_DATA(*,J))))

         INT_SIGNAL = INT_TABULATED(WAVELENGTH_DATA(X_INDEX_LEFT(0):X_INDEX_RIGHT(0),J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(X_INDEX_LEFT(0):X_INDEX_RIGHT(0),j))

         PRINT, 'INT_SIGNAL=' , INT_SIGNAL
         CGS_TITLE = 'INTEGRATED SIGNAL =' + STRING(INT_SIGNAL)
         LIMITS_CONTROL = 1

         ENDIF

         IF CONTROL EQ 3 THEN BEGIN
          spectra(*,j) = SMOOTH(spectra(*,j),3)
          PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), $
          XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
          TRANS_IMAGE_CUT(*,*,J) = SMOOTH(TRANS_IMAGE_CUT(*,*,J),3)
         ENDIF
         IF CONTROL EQ 4 THEN BEGIN
          spectra(*,j) = SMOOTH(spectra(*,j),5)
          PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), $
          XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
          TRANS_IMAGE_CUT(*,*,J) = SMOOTH(TRANS_IMAGE_CUT(*,*,J),5)
         ENDIF
         IF CONTROL EQ 5 THEN BEGIN
          spectra(*,j) = SMOOTH(spectra(*,j),7)
          PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), $
          XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
          TRANS_IMAGE_CUT(*,*,J) = SMOOTH(TRANS_IMAGE_CUT(*,*,J),7)
         ENDIF
         IF CONTROL EQ 6 THEN BEGIN
          spectra(*,j) = SMOOTH(spectra(*,j),9)
          PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), $
          XSTYLE=1, TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
          TRANS_IMAGE_CUT(*,*,J) = SMOOTH(TRANS_IMAGE_CUT(*,*,J),9)
         ENDIF


         IF CONTROL EQ 7 THEN BEGIN

          IF LIMITS_CONTROL EQ 0 THEN BEGIN

              MENU,['OK'],BLE,TITLE='MARK INTEGRAL LIMITS',BUTSIZE=400,xoffset=0,yoffset=0

          ENDIF ELSE BEGIN

              INT_SIGNAL_SLIT = FLTARR(512,NO_CAM)
              ;MENU,['1','2','3','4','5','7','8','9','10'],APERT_WIDTH,TITLE='CHOOSE APERTURE WIDTH IN PIXELS',BUTSIZE=400,xoffset=0,yoffset=0
              ARC_APERT = INDGEN(25)*0.2675 ; THIS IS VALUE OF ARCSEC PER 1 PXL ON SJ IMAGE

              MENU,['DONT PRESS THIS BUTTON',$
              ' 1pxl ='+STRING(ARC_APERT(1)),$
              ' 3pxl ='+STRING(ARC_APERT(3)),$
              ' 5pxl ='+STRING(ARC_APERT(5)),$
              ' 7pxl ='+STRING(ARC_APERT(7)),$
              ' 9pxl ='+STRING(ARC_APERT(9)),$
              '11pxl ='+STRING(ARC_APERT(11)),$
              '13pxl ='+STRING(ARC_APERT(13)),$
              '15pxl ='+STRING(ARC_APERT(15)),$
              '17pxl ='+STRING(ARC_APERT(17)),$
              '19pxl ='+STRING(ARC_APERT(19)),$
              '21pxl ='+STRING(ARC_APERT(21)),$
              '23pxl ='+STRING(ARC_APERT(23))],$
              APERT_WIDTH,TITLE='CHOOSE APERTURE WIDTH IN PIXELS',BUTSIZE=400,xoffset=0,yoffset=0

              APERT_WIDTH_recalc = APERT_WIDTH * 2 - 3

              FOR  PXL=0,511 DO BEGIN ;HERE J MEANS CAM INDEX
                 INT_SIGNAL_SLIT(PXL,J) = INT_TABULATED(WAVELENGTH_DATA(X_INDEX_LEFT(0):X_INDEX_RIGHT(0),J),$
                                        FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*TEMP_FOR_INTEGRATION(X_INDEX_LEFT(0):X_INDEX_RIGHT(0),PXL,j))

              ENDFOR

              ;IF APERT_WIDTH NE 1 THEN INT_SIGNAL_SLIT(*,J) = SMOOTH(INT_SIGNAL_SLIT(*,J),APERT_WIDTH)
              ;tutut

              IF APERT_WIDTH NE 1 THEN BEGIN
                 IF  APERT_WIDTH NE 2 THEN BEGIN
                   INT_SIGNAL_SLIT(*,J) = SMOOTH(INT_SIGNAL_SLIT(*,J),APERT_WIDTH_recalc)
                 ENDIF
              ENDIF ELSE BEGIN
                 C_D,1
              ENDELSE

              window,NO_CAM,ypos=80,xpos=80,xsize=320*4,ysize=128*5,title='VALUE OF INTEGRAL OF '+strupcase(CAM_INDEX(J))+' LINE PROFILE'
              PLOT, INT_SIGNAL_SLIT(*,J),TITLE = strupcase(CAM_INDEX(J))+' INTEGRATED SIGNAL'+DATA_TIME+' [UT]',$
              YTITLE='[ erg * sr^(-1) * cm^(-2) * s^(-1) ]',XTITLE='POSITION ON THE SLIT [ PXL ]',CHARSIZE=1.5,XSTYLE=1,$
              YSTYLE=1,YRANGE=[100.,MAX(INT_SIGNAL_SLIT(*,J))], /YLOG

          ENDELSE

         ENDIF

         IF CONTROL EQ 8 THEN BEGIN
          REPEAT BEGIN
              MENU,['OK','+100','+10','+1','+0.1','-0.1','-1','-10','-100'],CORR_INT,TITLE='OFFSET CORRECTION',BUTSIZE=300,xoffset=0,yoffset=-15
              ;VECTOR OF VALUES;
              CORR_INT=CORR_INT-1
              FAC_INT= [0.,100.,10.0,1.,0.1,-0.1,-1.0,-10.0,-100.]
              spectra(*,j) = spectra(*,j) + FAC_INT (CORR_INT)
              TEMP_FOR_INTEGRATION(*,*,J) = TEMP_FOR_INTEGRATION(*,*,J) + FAC_INT (CORR_INT)
              PLOT,WAVELENGTH_DATA(*,J),FACTOR_CGS(J)*FACTOR_DATA_VECTOR(*,J)*spectra(*,j), XSTYLE=1, YSTYLE=1,$
              TITLE=strupcase(CAM_INDEX(J))+' SPECTRUM ',charsize=1.5,YTITLE=Y_TITLE(J)
              PLOTS,[MIN_WL(J),MAX_WL(J)],[0,0],/DATA,COLOR=150

          ENDREP UNTIL CORR_INT EQ 0


         ENDIF
         IF CONTROL EQ 9 THEN BEGIN
          MENU,['ASCII', 'FITS'],SAVE_CONTROL,TITLE='FILE FORMAT',BUTSIZE=300,xoffset=0,yoffset=-15


          IF SAVE_CONTROL EQ 1 THEN BEGIN
              OPENW,2,'flux_integral_'+data_time(0)+'_'+strlowcase(CAM_INDEX(J))+'.txt'
              FOR PXL=0,511 DO BEGIN
                 PRINTF,2,X_SLIT_POSITION,PXL,INT_SIGNAL_SLIT(PXL,J)
              ENDFOR
              CLOSE,2
          ENDIF
          IF SAVE_CONTROL EQ 2 THEN BEGIN
              ;WRITEFITS INT_SIGNAL_SLIT(*,J)

              WRITE_MATRIX = FLTARR(3,512)
              WRITE_MATRIX(0,*) = X_SLIT_POSITION
              WRITE_MATRIX(1,*) = PXL
              WRITE_MATRIX(2,*) = INT_SIGNAL_SLIT(*,J)

              writefits, 'flux_integral_'+data_time(0)+'_'+strlowcase(CAM_INDEX(J))+'.fit', WRITE_MATRIX
          ENDIF
         ENDIF

       ENDREP UNTIL CONTROL EQ 2


    ENDIF

    IF (LETR EQ 15) THEN BEGIN  ;SHOW OTHER IMAGES
       ;MENU,DATA_TIME,FIRST_FILE_INDEX,TITLE='CHOOSE  .:: THE FIRST ::.  FILE ',BUTSIZE=400
        ;IMAGE_INDEX = FIRST_FILE_INDEX-1    ;CHANGING GLOBAL VALUE
        ;MENU,DATA_TIME(FIRST_FILE_INDEX-1:*),LAST_FILE_INDEX,TITLE='CHOOSE  .:: THE LAST ::.  FILE ',BUTSIZE=400

        FIRST_data_name = dialog_pickfile(title='CHOOSE THE FIRST FILE', /file,  /MULTIPLE_FILES, get_path = dir, FILTER = ['*sj*_R.fit'])
       filelist_SHOW = FINDFILE(dir+'*sj*_R.fit')
       FIRST_FILE_INDEX = WHERE (FIRST_data_name(0) EQ filelist_SHOW)

        LAST_data_name = dialog_pickfile(title='CHOOSE THE LAST FILE', /file,  /MULTIPLE_FILES, path = dir, FILTER = ['*sj*_R.fit'])
       filelist_SHOW = FINDFILE(dir+'*sj*_R.fit')
       LAST_FILE_INDEX = WHERE (LAST_data_name(0) EQ filelist_SHOW)

       II=0
       FOR JJ=FIRST_FILE_INDEX(0) ,LAST_FILE_INDEX(0) DO BEGIN

         SHOW_DATA = readfits(filelist_SHOW(JJ),header_SHOW)
         window,II,ypos=(II/4)*150,xpos=(II MOD 4) *325,xsize=320,ysize=128,title=strcompress(HEADER_SHOW(6),/remove_all)
         tvscl,REBIN (SHOW_DATA,320, 128)
         II=II+1
         WAIT, 0.06
       ENDFOR
       MENU,['OK'],C,TITLE='SJ IMAGES',BUTSIZE=400
       FOR U=0,II DO WDELETE,U & WAIT,0.1

    ENDIF


    IF (LETR EQ 16) THEN BEGIN   ;SAVING SPECTRA

        MENU,['ASCII', 'FITS'],SAVE_CONTROL,TITLE='FILE FORMAT',BUTSIZE=300,xoffset=0,yoffset=-15

          IF SAVE_CONTROL EQ 1 THEN BEGIN
              FOR JJ=1,NO_CAM-1 DO BEGIN
                 OPENW,2,'spectrum_'+data_time(0)+'_'+strlowcase(CAM_INDEX(JJ))+'.txt'
                 FOR PXL=0,1279 DO BEGIN
                   PRINTF,2,WAVELENGTH_DATA(PXL,JJ),FACTOR_EXP(JJ)*FACTOR_CGS(JJ)*FACTOR_DATA_VECTOR(PXL,JJ)*spectra(PXL,JJ)
                 ENDFOR
                 CLOSE,2
              ENDFOR
          ENDIF
          IF SAVE_CONTROL EQ 2 THEN BEGIN
              FOR JJ=1,NO_CAM-1 DO BEGIN
                 WRITE_MATRIX = FLTARR(2,1280)
                 WRITE_MATRIX(0,*) = WAVELENGTH_DATA(*,JJ)
                 WRITE_MATRIX(1,*) = FACTOR_EXP(JJ)*FACTOR_CGS(JJ)*FACTOR_DATA_VECTOR(*,JJ)*spectra(*,JJ)
                 writefits, 'spectrum_'+data_time(0)+'_'+strlowcase(CAM_INDEX(JJ))+'.fit', WRITE_MATRIX
              ENDFOR
          ENDIF
    ENDIF


ENDREP UNTIL LETR EQ 3


for i=0,no_cam do wdelete,i



print,'END OF JOB'

end