AWARES ExOb Wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Jans IDL tips ======

Dies ist eine Sammlung meiner idl //Entdeckungen//
This is a collection of my IDL //discoveries//
\\

In the web you may find many IDL tutorials or helb pages ...\\

IDL current reference @ NV5 Geospatial software :\\
[[https://www.nv5geospatialsoftware.com/docs/routines-1.html| IDL  8.9. Reference @ NV5 ]]\\
[[https://www.nv5geospatialsoftware.com/docs/funclisting.html| IDL  8.9. funclisting @ NV5]]\\
[[https://www.nv5geospatialsoftware.com/docs/WhatsNew.html| IDL  8.9. whats new @ NV5]]\\



[[https://www.l3harrisgeospatial.com/docs/routines-1.html | IDL  8.9. references @ l3harris geospatial ]]\\

[[http://www.harrisgeospatial.com/docs/routines-1.html | IDL 8.8 referecnes ]]  (dead, Certifacte expired ??? )
It was once Exelis, but they deviate to harris:\\
[[http://www.exelisvis.com/docs/routines-1.html  |  refefrences ]] ... links to newest version at harris \\

A site which provides the IDL-help pages which you also cen get at the IDL prompt with %%IDL>? %%//functionname//%% %%\\
[[https://www.irya.unam.mx/computo/sites/manuales/IDL/idl_Left.html]] 


At FZ-Jülich you can find a manual for IDL programming:\\
[[http://www2.fz-juelich.de/jsc/files/docs/vortraege/idl.pdf|Praxisbezogene IDL Programmierung]]


You can get Help from IDL by just typing ? at the IDL prompt
or for a certain idl function by typing %%? <functioname>%%.
This will open your web browser with a local help page.


We once decided to use a **[[internal:administration:rcs|revision control system (RCS)]]** for the idl routines under /home/hatpro/idl


===== Configuration =====

IDL uses some Environment variables containing the pathes to its libraries etc. see [[https://www.nv5geospatialsoftware.com/docs/prefs_directory.html | here ]].
Internally some of them can be accessed as system variables: The environment varibles usually have the form IDL_<NAME>. Inside IDL they are accesible via !<NAME>
Basic path to IDL is in environamet variable ''IDLDIR'' inside IDL it is accesible via ''!DIR'' which is currently (idl version 9.0) ''/opt/nv5/idl90'' .



===== GDL =====

There is an alternative for IDL based on the same syntax, function names etc.: \\
Gnu Data Language (**GDL**)\\

[[http://gnudatalanguage.sourceforge.net/]]\\
[[http://linux.die.net/man/1/gdl]]\\
[[http://blog.redbranch.net/2011/04/13/installing-gdl-an-idl-alternative/|weiterführende hilfe für die installation]]\\
[[http://aramis.obspm.fr/~coulais/IDL_et_GDL/Matrice_IDLvsGDL_missing.html|list of functions missing in GDL]]\\
... auf der qeb ist das standard ubuntu paket //gnudatalanguage// version 0.9 installiert.\\
Call with\\ 
%%$>gdl%%\\
But it is since years in development and many functions are missing or not working properly.\\
See this list: %%/home/jschween/ideen/gdl/readme.txt%%\\

Starten eines Scirptes von der Shell geht in GDL 0.9 nicht wie in IDL mit %% \$> idl @idl_script %% sondern über eine bash pipe:\\ %% \$> cat @idl_script | gdl %%

Einige Besonderheiten:
siehe diese [[http://www.sgeier.net/tools/GDL-intro.pdf | einfachste Einleitung in GDL]]
  * shell Befehle müssen nicht per %%spawn%% gestartet werden sondern können in einer Zeile die mit %%\$%% beginnt stehen, z.B.:\\ %%GDL> $ cp idl_script.pro gdl_script.pro%%
  * ...






===== Basics =====

start from the terminal prompt:

%%  ~>idl%%

Help with

%%  IDL> ?%%

or

%%  IDL> ? command-name%%


idls own routines are located in 

%%  /usr/local/itt/idl64/lib%%


==== simple examples ====

=== hello world ====

%%  print, 'hello world !'%%


=== plot ====

<code idl>
  ; number of data points
  N = 100
  ; an array of x-values in interval 0..1
  x = findgen(N)/N
  ; an array of y-values
  y = sin(x*2*!Pi)
  ; plot the data
  plot, x, y
</code>

add axis titles and title to plot:

<code idl>
  plot, x,y , title='this is the title' , xtitle='x-axis', ytitle='ordinate'
</code>

manipulate the range

<code idl>
  plot, x,y , xrange = [ -0.1 , 1.5 ] , xstyle=1 , yrange = [ -1.2 , +1.5 ] , ystyle=1
</code>

use linestyle and colors

<code idl>
  plot, x,y , linestyle=1, color=3
</code>

plot a second function over the first

<code idl>
  oplot, x, 0.6*sin(2*x) , linestyle=2, color=255
</code>


=== export data to file ====

<code idl>
  openw, file, 'data.txt',  /Get_LUN
  printf, file , 'x=', x 
  printf, file , 'y=', y
  close, file
  free_lun, file
</code>

=== calculus ====

<code idl>
  ; add two to all values
  x2 = x+2
  ; multiply all by three
  x3 = x*3
  ; square and squareroot
  x_square = x^2
  x_root = sqrt(x)
  ; mean etc.
  ym = mean(y)
  yvar = variance(y)
  ysigma = stddev(y)
 </code>


==== general syntax: ====

  <command> [ , <arguments> ] [ , <name>=<variable> ] [ , /<flagname> ]

  * be aware of the comma !
  * lowercase/uppercase does not matter
  * <arguments> is a comma separated list like x,y,z which must come in the specified order
  * the last arguments can be omitted, e.g. x,y instead of x,y,z (if <command> can handle this )
  * <name>=<variable> are named variables inside command, these parameters can usually be omitted
  * <name>  arguments must not come in a special order
  * /<flagname> is an abbreviation for flagname=1


==== change working dir ====

  cd , "<path>" [ ,current=<var> ]

path in single (') or double (") quotation marks !

==== execute shell command : SPAWN ====

  SPAWN [, "<Command>" [, Result] [, ErrResult] ]

command in single (') or double (") quotation marks !

e.g. pwd = print working directory:

  IDL>SPAWN , 'pwd'

In some cases you may get error message from your program. For example curl:
  IDL> spawn, 'curl -V'
  curl: /opt/nv5/idl90/bin/bin.linux.x86_64/libcurl.so.4: no version information available (required by curl)
  ...
  WARNING: curl and libcurl versions do not match. Functionality may be affected.

The first and the last lines tell you that %%libcurl.so.4%% from the %%/opt/nv5/idl90/...%% path makes problems. 
But why should curl try to load libraries from /opt/nv5/... i.e. the IDL folder ?

Another example is ghostscript (gs) - which does not even work:
  IDL> SPAWN, 'gs --version'
  gs: /opt/nv5/idl90/bin/bin.linux.x86_64/libstdc++.so.6: version `GLIBCXX_3.4.32' not found (required by /lib/x86_64-linux-gnu/libgs.so.10)

Again ghostscript tries to load lubraries from /opt/nv5/... 

To get around you have to understand where the system seeks libraries.
This is explained e.g. here: [[https://askubuntu.com/questions/1500913/find-default-library-paths-location-in-ubuntu]], or in the man pages of [[https://man7.org/linux/man-pages/man8/ld.so.8.html | ld.so]]. To ensure that IDL only uses libraries of a certain version they set the variable LD_LIBRARY_PATH to the idl library folders. It is not set in standard bash environment, but in the SPAWN environment as it is a subrpocess of IDL:
  IDL> SPAWN, 'echo $LD_LIBRARY_PATH'
  /opt/nv5/idl90/bin/bin.linux.x86_64:/opt/nv5/idl90/bin/bin.linux.x86_64/dm/lib:/opt/nv5/idl90/bin/bin.linux.x86_64/jre/lib/server:/opt/nv5/idl90/bin/bin.linux.x86_64/jre/lib/server/..:/opt/nv5/idl90/bin/bin.linux.x86_64/jre/lib/server/../native_threads

We can ask Linux where we can find the libraries needed:
  $whereis libcurl.so.4
  libcurl.so.4: /usr/lib/x86_64-linux-gnu/libcurl.so.4
and
  $whereis libstdc++.so.6
  libstdc++.so.6: /usr/lib/x86_64-linux-gnu/libstdc++.so.6

We have several possibilities:

  * adapt the LD_LIBRARY_PATH from within IDL before invoking SPAWN. This was proposed to me : [[https://www.nv5geospatialsoftware.com/Support/Forums/aft/8650 in the IDL forum]] for the PAtH variable but maybe IDL will have problems with that. And you should do this only once to avoid an extremly long path with only your prepended entry.
      %%SETENV, 'LD_LIBRARY_PATH=/usr/lib/x86_64-linux-gnu/:'+GETENV('LD_LIBRARY_PATH')%%

  * Alternatively adapt LD_LIBRARY_PATH within the SPAWN environment when calling your program
      %%SPAWN, 'LD_LIBRARY_PATH=/usr/lib:$LD_LIBRARY_PATH; your_command_here'%%

  * or even simpler unset this path variable in SPAWN: it is not needed as we are not doing IDL specific things, and when SPAWN is finished this is forgotten again.
      %%SPAWN, 'unset LD_LIBRARY_PATH; your_command_here'%%

According to IDL help about [[https://www.nv5geospatialsoftware.com/docs/Environment_Variables.html | Environment variables]] LD_LIBRARY_PATH is used for the Java libraries - i.e. probably for the IDE and its graphical interface. That means if you do not use the IDE the first solution will have no problems.

Check with curl:
    SPAWN, 'unset LD_LIBRARY_PATH ; curl -V'
    curl 8.5.0 (x86_64-pc-linux-gnu) libcurl/8.5.0 OpenSSL/3.0.13 zlib/1.3 brotli/1.1.0 zstd/1.5.5 libidn2/2.3.7 libpsl/0.21.2 (+libidn2/2.3.7) libssh/0.10.6/openssl/zlib nghttp2/1.59.0 librtmp/2.3 OpenLDAP/2.6.7
    ...

Check with ghostscript
    SPAWN, 'unset LD_LIBRARY_PATH ; gs --version'
    10.02.1 

=> WORKS!

==== compile and run idl programs ====

to just compile do

  .COMPILE <filename_1> [,<filename_2> [...]]

to compile and run:

  .RUN <filename_1> 

with <filename_1> without the extension ".pro"


==== general structure ====

The semicolon ';' is the comment sign: everything after it is ignored by the interpreter

a PROGRAM is a text file with IDL commands ending with an END statement (lower and  uppercase is allowed, dont mix cases !)
  print, 'Hello world !'
  print, '3*4=', 3*4
  x = 3
  y = 2.5
  z = x/y
  print, 'x=', x, ' y=',y,' x/y=z=', z
  print, 'its time to say goodbye'
  END


decisions are made with IF THEN ELSE 

  IF x eq 4 THEN print, 'x equals 4' ELSE print, 'x is not 4'
  
this must be in one line but can be extended of several lines with the '$' character:
  
  IF x eq 4 $
    THEN print, 'x equals 4' $
    ELSE print, 'x is not 4'

or with BEGIN END blocks

  IF x eq 4 THEN BEGIN
    print, 'x equals 4' $
  END ELSE BEGIN
    print, 'x is not 4'
  END

for safety you may extend the END statements with blocked statement i.e. ENDIF or ENDELSE

you may nest i.e. build structures like IF ... THEN ... ELSE IF ... THEN ...
or use the CASE statement.

Loops can be implented with FOR, WHILE or REPEAT as one-liners:

  FOR i = 0,5 DO PRINT, i

  i = 0  
  WHILE i lt 6 DO i = i+1
  
  i = 1
  REPEAT i = i + 1 UNTIL i gt 5

for more complex loops you may use BEGIN ... END

  FOR i = 0,5 DO BEGIN
    PRINT, i
    END
  
  i = 0
  WHILE i lt 6 DO BEGIN
    i = i+1
    PRINT, i
    END
  
  i = 1
  REPEAT BEGIN
    i = i + 1 
    PRINT, i
  END UNTIL i gt 5

to make coding a bit safer you may append the name of the loop command to the END statement, i.e.
  ENDFOR 
  ENDWHILE
  ENDREP 

a PROCEDURE is a subroutine that may take parameters and change them. Its definition is started with the word PRO and ends with a END statement

  PRO my_proc, text_str, x_value=x_value
    print, 'text_str=', text_str
    print, 'x=', x_value
    x_value = x_value + 2
    END ; pro

it is called with its parameter list separated by commas:
  x = 2
  my_proc, 'some text', x_value=x
  
In this case x will have the value x+2 after the call 

if a procedure is stored in a file with the same name as the procedure IDL will compile this file ...


A FUNCTION is a subroutine that may return a result

  FUNCTION my_func, text_str, x_value=x_value
    print, 'text_str=', text_str
    print, 'x=', x_value
    x_new = x_value + 2
    RETURN, x_new
    print, 'this is never printed !'
    END ; function

return means that the result is returned but also that the execution returns from the sub subroutine, for that reason the last print command is ignored.

A function is called with its parameterlist in braces

  x = 2
  y = my_func( text_str, x_value=x )
  

==== idl path ====

When finding a new name which might be a function procedure or another idl program idl searches in the directories given by the idl sytem variable **!PATH**. To make own programs visible to idl the directory where these programs could be found must be **included in the path**:

  !PATH = !PATH + ':/home/hatpro/idl/lib'

To ensure that the directory is not added again an again use the following code:

  s = '/home/hatpro/idl/lib'
  if strpos( !PATH , s ) lt 0 then !PATH = !PATH + ':' + s


==== defining constants ====

Constants are defined by value, the type can be specified
  <name> = <Value>[<Type Letter[s]>]

e.g.:

  * i = 0b ;     creates a      1byte=8bit integer constant of value zero (0..255) 
  * i = not 0b ; creates a            8bit integer with all bits set to 1 i.e. value 255decimal='ff'x='11111111'b
  * i = 0S ;     creates a    signed 16bit integer constant of value zero (-32767..32767) 
  * i = 0US ;    creates an unsigned 16bit integer constant of value zero (0..65535)
  * z = 0L ;     creates a    signed 32bit integer constant (Long) of value zero (+/-2.2mrd)
  * z = 0UL ;    creates an unsigned 32bit integer constant (Long) of value zero (0...4.mrd)
  * z = 0LL ;    creates a    signed 64bit integer constant (Long) of value zero (+/-9.2e18)
  * z = 0ULL ;   creates an unsigned 64bit integer constant (Long) of value zero  (0...1.8e19)
  * x = '0F'XU ; an unsigned integer (U) written as Hex (X) with (decimal) value 15
  * x = '08'OU ; an unsigned integer (U) written as Octal (O) with (decimal) value 8
  * y = 7.3E1 ; creates single precision 73.0
  * z = 7.3D1 ; creates double precision 73.0
  * novalue = !VALUES.F_NAN ; creates a single precision floating point NAN ('not a number')
  * novalue = !VALUES.D_NAN ; creates a double precision floating point NAN ('not a number')
  * infty = !VALUES.F_INFINITY ; creates a single precission infinity (positive, a minus in front makes it negative)
  * infty = !VALUES.D_INFINITY ; creates a double precission infinity
...

see [[https://www.harrisgeospatial.com/docs/Defining_and_Using_Const.html|'defining constants' ]] and [[https://www.harrisgeospatial.com/docs/IDL_Data_Types.html|'data types]] in IDL help. 

For NAN and infinity see 'Special Floating-Point Values' or [[https://www.harrisgeospatial.com/docs/Constant_System_Variable.html|Constant System Variables]]

single precission is a 4 Bit floting point number (~range: -1e38..1E38 accuracy ~1e-7),\\
double precision is a 8Bit floating point number (~range: -1E308..1E+308 ) \\
see [[http://en.wikipedia.org/wiki/Floating-point|wikipedia floating point]].

==== arrays ====

can be generated with 
  x1d = fltarr( N ) ; type float
  x2d = dblarr( Ni, Nj ) ; type double precision
  ...
The arrays will be filled with zeros, see e.g. help for [[http://www.physics.nyu.edu/grierlab/idl_html_help/F29.html#wp676939| fltarr ]] or the info page about [[http://www.physics.nyu.edu/grierlab/idl_html_help/funclisting4.html#wp1154150|array creation ]]

To fill the array with a certain value:
  x2d[*] = -1.0

To fill it with the value of its running index:
  x1d = findgen( N ) ; float type
  x2d = sindgen( Ni, Nj ) ; string type
  ...
For details see [[http://www.physics.nyu.edu/grierlab/idl_html_help/F22.html#wp676873|here]]

You may also use [[http://www.physics.nyu.edu/grierlab/idl_html_help/M3.html#wp887420|makearray]]:
  s1d = MAKE_ARRAY( Ni, /string, /index )
  x2d = MAKE_ARRAY( Ni, Nj , /float, value=-1.0 )
  ...
  
Multidimensional arrays can be accessed via indices for all the different dimensions but also via less indices:

  x3d = fltarr( Ni, Nj, Nk )
  print, x3d[i1d] ; returns the i1d scalar element of the array
  print, x3d[i,j] ; returns the i+j*Ni scalar element of the array
  print, x3d[i,j,k] ; returns the i+j*Ni+k*Ni*Nj((correct ?)) scalar element of the array
  ...

For two dimensional arrays the relation between one and two dimensional indices is:

  i1d = i + j*Ni
  i = i1d mod Ni
  j = i1d / Ni
  
This can also be done with the idl function [[https://www.l3harrisgeospatial.com/docs/ARRAY_INDICES.html|ARRAY_INDICES]]


==== structures ====

A structure is a dataype which joins variables of different type in one variable.
It can be generated with (see [[http://www.harrisgeospatial.com/docs/Creating_and_Defining_St.html]])

  var = { name : 'bill' , int1:1 , a2:2.4 }
  
or with the CREATE_STRUCT function (see [[http://www.harrisgeospatial.com/docs/create_struct.html]])

  var = CREATE_STRUCT( 'name', 'bill' , 'int1' , 1 , 'a2' , 2.4 }

usage is
  var = CREATE_STRUCT(  [Name=<string>, ] <name1>, value1, [<name2>, value2, ... ] ...    )
where Name registers this strucure to IDL and forces you to not deviate from the list of names and types ...
or
  var = CREATE_STRUCT(  [Name=<string>, ] <namelist>, value1, [ value2, ... ] ...    )
with <namelist> a string array e.g. <namelist> = [ 'name', 'int1', 'a2' ] 

CREATE_STRUCT allows to append structures
  var = CREATE_STRUCT(  [Name=<string>, ] [struct_1, ...] <name1>, value1, [<name2>, value2, ... ]  [struct_2,...]    )
Where struct_1, and struct_2 are previosly defined structures.

reference is made with 

  print, var

or with the . operator:

  print, var.name
  print, var.int1
  ...

if you do not know the elements of a strucuture you may use help with the flag /strucutre

  help, var, /structure

or the TAG_NAMES function

  tn = TAG_NAMES(var)
  print, tn
  
You can construct arrays of structures with the REPLICATE function and manipulate its content and find certain elements (see [[http://www.harrisgeospatial.com/docs/Arrays_of_Structures.html]]):

  ; create array of 20 structures var as e.g. defined above:
  var_list = REPLICATE( var , 20 )
  ; change name of entry name of structure number three
  var_list[3].name='ff'
  ; find this entry
  print, where(var_list.name eq 'ff' )
     

==== LIST and HASH ====

With version 8.0 IDL introduced LIST and HASH.
  * LIST is an array like datatype which may contain a mixture of different variable types and structures. Elements can be accessed via their index (see [[https://www.irya.unam.mx/computo/sites/manuales/IDL/Content/Reference%20Material/L/LIST.html|LIST help page at irya.unam.mx]])
  * HASH is a structure like dataype which elements can be accessed via strings containing the names and which is more flexible than the STRUCUTRE data type (see [[https://www.irya.unam.mx/computo/sites/manuales/IDL/Content/Reference%20Material/H/HASH.html|HASH help page at irya.unam.mx]])

In the following %%value%% can be of any IDL datatype.


Create a LIST:
  x_list = LIST( Value1, Value2, ... )
Add a single element
  x_list.ADD, value
  ...
  
Create a HASH:
  y_hash = HASH( [Key1, Value1, Key2, Value2, ... )
Add another HASH
  y_hash = y_hash + HASH( [Key21, Value21, ... )
  ...
  
Both datatypes are objects with several methods to get information or manipulate its contents - se help page referenced above.

The FOREACH statement can be applied to both data types. Assuming that you can print all their elements you can do this:
  FOREACH x in x_list do print, x
  FOREACH y in y_hash do print, y
  ...
  

==== matrix multiplication ====

IDL comes with two matrix multiplication operators %%#%% and %%##%%.
  * lets assume we have A=fltarr(NAcol,NArow) and B=fltarr(NBcol,NBrow) 
  * operator %%#%% multiplicates rows of the left by columns of the right matrix \\ => it must be NArow=NBcol and we get C = A # B = fltarr(NAcol,NBrow) \\ it is %%C[i,j] = sum(k=0..NArow-1) A[i,k] * B[k,j]%% \\ scheme:

                            bbb
                  aaaaa     bbb     ccc
  C  =  A # B  =  aaaaa  #  bbb  =  ccc
                  aaaaa     bbb     ccc
                            bbb

  * operator **%%##%%** multiplicates columns of the left by rows of the right matrix \\ => it must be NAcol=NBrow and we get C = A ## B = fltarr(NBcol,NArow) \\ it is C[i,j] = sum(k=0..NAcol-1) A[k,i]*B[j,k] \\ scheme:

                   aaa               cccc
                   aaa      bbbb     cccc
  C  =  A ## B  =  aaa  ##  bbbb  =  cccc
                   aaa      bbbb     cccc
                   aaa               cccc 

  * you may use this to calculate the scalar product of two vectors: \\ s = transpose(x) # y  or  s = x ## transpose(y)  with x = fltarr(N) and y = fltarr(N) 

  * IDL is somewhat sloppy with the relation between vectors and matrices, i.e. one and two dimensional arrays:
    * lets assume X = fltarr(N) and Y = fltarr(N) 
    * the following gives the same N*N matrix:
      * C = X # Y
      * C = X # transpose(Y)
    * and these give the transpose of C (D = transpose(C))
      * D = X ## Y
      * D = X ## transpose(Y)

  * you may generate Matrices from vectors, i.e. from one-dimensional arrays:
  * with X=fltarr(Nx) and Y = fltarr(Ny) you get (again using the flopiness of IDL)
    * C = X # Y  or  C = X # transpose(Y)   => C = fltarr(Nx,Ny)
    * C = X ## Y  or C = transpose(X) ## Y => C = fltarr(Ny,Nx)
  * in linear algebra a product between two vectors giving a matrix is called the [[https://en.wikipedia.org/wiki/Outer_product|'outer product']] or in german [[https://de.wikipedia.org/wiki/Dyadisches_Produkt|'Dyade']]

  * you may use this to generate matrices with equal rows or columns:
    * to generate a Matrix with Ny equal rows    x = fltarr(Nx) use \\ XX = x # replicate(1,Ny) \\ i.e. replicate(1,Ny) is interpreted as a 1*Ny matrix
    * to generate a Matrix with Nx equal columns y = fltarr(Ny) use \\ YY = replicate(1,Nx) # y \\ i.e. replicate(1,Nx) is interpreted as a Nx*1 matrix - odd, he ! - try to be not confused !

                           1     xxx
                           1     xxx
  XX  =  X # E  =  xxx  #  1  =  xxx
                           1     xxx
                           1     xxx


                           y     yyy
                           y     yyy
  YY  =  E # Y  =  111  #  y  =  yyy
                           y     yyy
                           y     yyy 

  * One can use this to generate matrices with a variable increasing only in one dimension and arrays of function values:

  ; generate Nx*Ny Matrix xx with Nx values increasing in x-direction from xa to xb 
  x = xa + (xb-xa) * findgen(Nx)/(Nx-1) 
  xx = x # replicate(1.0,Ny)
  
  ; generate Nx*Ny Matrix yy with values increasing in y-direction from ya to yb
  y = ya + (yb-ya) * findgen(Ny)/(Ny-1) 
  yy = replicate(1.0,Nx) # y
  
  ; Nx*Ny matrix with distances from a point with coordinates (cx,cy)
  r = sqrt( (xx-cx)^2 + (yy-cy)^2 )
  
  ; Nx*Ny matrix with function values on every point 
  f_xy = f( xx, yy )

  * As matrix operations are faster than for-loops this may speed up your code a lot ...

==== procedure names ====

if an unknown procedure name <fname> is found idl looks in 
  - the current directory
  - in the directories specified in IDLs '!PATH'-variable

for a file <fname>.pro or  <fname>.sav

letters of the filename should be all lowercase.

the !PATH variable ist constructed as 
  <dir1>:<dir2>:...


==== function names ====

if idl encounters an expression like
  y = f_phi2( a )

It assumes firstly an array variable f_phi2 of which the index a is desired.\\
If f_phi2 is an external function it should be anounced with the command 
  FORWARD_FUNCTION f_phi2

see also  
http://www.physics.nyu.edu/grierlab/idl_html_help/F33.html#wp872909



==== Blocks ====
 
Blocks are generated with:
  BEGIN
  ...
  END
where END can expanded to
  ENDIF   	
if it terminates an IF THEN ... block
  ENDWHILE
if it terminates a WHILE DO ... block
...

This naming of the END key word is helpful to find errors in the BEGIN...END structure 
==== Interaktion ====

 
soll laut idl mit //widgets// gemacht werden - ist aber aufwändig !!

=== Tastatur ===

es gibt zwei einfache methoden mit der Tastatur zu Arbeiten:

  * der **Read**-Befehl erwartet eine Eingabe die mit RETURN abgeschlossen werden muss
  * der **GET_KBRD** erlaubt eine direkte Tastaturabfrage:

der **Read**-Befehl:
<code idl>
c = ' '
read,'continue with RETURN',c
</code>

mit **GET_KBRD** kann man Tastaturabfragen machen:
<code idl>
  Result = GET_KBRD([Wait] [, /ESCAPE | , /KEY_NAME])
</code>

Zum Beispiel:
  key = GET_KBRD( 0 )
liefert in key die letzte gedrückte Taste oder %%''%% (nullstring) wenn keine Taste gedrücktwurde ...

  key = GET_KBRD( 1 )
wartet auf tastendruck und liefert in key die letzte gedrückte Taste

  key = GET_KBRD( 1 , /escape )
wartet auf tastendruck und liefert in key die letzte gedrückte Taste oder 
die escape sequenz die der Taste entspricht - also zb. string(byte([27,91,66])) 
für die 'pfeil nach unten' taste

Liste verschiedener escape sequenzen:
|  Taste             |  escape sequence |\\
| cursor n. oben     | 27 91 65 |
| cursor n. unten    | 27 91 66 |
| cursor n. rechts   | 27 91 67 |
| cursor n. links    | 27 91 68 |
| ende               | 27 91 70 |
| pos 1              | 27 91 72 |
| Einfg              | 27  91  50 126 |
| Entf               | 27  91  51 126 |
| seite n. oben      | 27  91  53 126 |
| seite n. unten     | 27  91  54 126 |
| ESC                | 27  27 |


damit lässt sich ein einfache block aufbauen mit dem man ein Programm steuern kann:\\
  * das Programm soll eine variable i hoch zählen
  * anhalten wenn %%' '%% (space) gedrückt wird
  * eins weiter (i+1) gehen wenn '>' oder 'cursor nach rechts' gedrückt wird
  * eins zurück (i-1) gehen wenn '<' oder 'cursor nach links' gedrückt wird
  * an den anfang (i=0) gehen wenn '0' oder 'pos 1' gedrückt wird
  * an das Ende (i=N-1) gehen wenn 'E' oder 'Ende' gedrückt wird
  * beenden wenn ESC gedrückt wird

<code idl>
i = 0
finished = 0
key = ' '
while not finished do begin

  ; ... do something ...
  print, i,' ... '

  ; check keyboard - return escape sequences 
  ; if the last key was a space - wait
  if key ne ' ' then key = GET_KBRD( 0 , /ESCAPE ) $
                else key = GET_KBRD( 1 , /ESCAPE )

  ; if we got a key
  if key ne '' then begin
    if key eq ' ' then begin
      ; ' '(space) => increment counter
      if i lt N-1 then i += 1
    endif else if (key eq '<') or (key eq string(byte([27, 91, 68]))) then begin
      ; '<' or cursor left => go back
      if i gt 0 then i -= 1
      key = ' '
    endif else if (key eq '>') or (key eq string(byte([27, 91, 67]))) then begin
      ; '>' or cursor right => go further
      if i lt N-1 then i += 1
      key = ' ' 
    endif else if (key eq '0') or (key eq string(byte([27, 91, 72]))) then begin
      ; '0' or pos 1 => go to start
      i = 0
      key = ' ' 
    endif else if (key eq 'e') or (key eq string(byte([27, 91, 70]))) then begin
      ; 'e' or Ende => go to end
      i = N-1
      key = ' ' 
    endif else if key eq string(byte([27,27])) then begin
      ; ESC-key => set stop flag
      finished = 1
    endif
  endif else begin ; no key pressed
    ; increment counter
    if i lt N-1 then i += 1
    ; set 'stop flag' if at end of data
    if i eq N-1 then key = ' '
    print,'End of data ...'
  endelse

endwhile 
</code>

Dieser code findet sich als funktion %%key_control.pro%% in %%hatpro/idl/lib%%.

=== mouse ===

mit der funktion **%%cursor%%** es ist möglich die Mausposition abzufragen:
<code idl>
  ; wait for mousebutton down event
  cursor, mx,my, /down
</code>
oder
<code idl>
  ; get mouse position imediately - point coordinates in data-coordinates
  cursor, mx,my, 0, /data
</code>

Information über die gedrückte Maustaste erhält man über die systemvariable %%!Mouse.button%% (linke maustaste == 1, rechte ==2, mitte == 4).

Eine Anwendung findet sich in **%%mouse_clip.pro%%** in %%hatpro/idl/lib/%%.

==== strings ====

can be split by
  s2 = STRMID( s , index , [length , ... ] )

position of substrings with 
  i = STRPOS( s , sub [ ... ] )

split string at character sep into substrings, return array of substrings (sx) and give number of substrings (N)
  sx = STRSPLIT( s , sep , count=N, /extract )
can be used to separate columns - see text files below...

you can also use regular expressions - e.g. lets assume you have strings with variable names and unit in braces:
  s = 'wind speed (m/s)'
then you may want to split at the opening and closing braces:
  ss = strsplit( s, '(\(|\))', /extract, /regex )
will deliver
  ss[0] = 'wind speed'
  ss[1] = 'm/s'


join array of strings into one string with separating character sep
  s = STRJOIN( s , sep )


filenames can be separated into path and name by
  file_dirname( <path to file> )
  file_basename( <path to file> )

==== formatted output ====

Formatted output is generated with a syntax similar to FORTRAN:

  PRINT, 'ABC..', 1, 2, 1.2, FORMAT='(A8,I," / ",I03," * ",F5.1)'
generates:
   ABC..      1  / 002 *   1.2
see widhts and positions:
   123456781234  / 123 * 12345

where the %%FORMAT='( ... )'%%  specifier said:
  * A8 : print as alpha-numerc-string of 8 characters width
  * I  : print as integer with standard width (8 digits, leading blanks)
  * " / ": insert string ' / '
  * I03: print as integer with 3 digits width and leading zeroes
  * " * ": insert string ' * '
  * F5.1: print as flaoting point with total width 5 and one digit after the decimal point

To print an array x=fltarr(N) with N ≥ 10 into a table with colums separated by tabs (ascii code 9) use:

  TAB=STRING(9b)
  PRINT, 'x=',x,FORMAT='(A,10("'+TAB+'",F8.2))'

Time as Julian date can be formatted as well. To get a date as DD.MM.YYYY hh:mm:ss use:

  PRINT, t, FORMAT='(C(CDI02,".",CMoI02,".",CYI04," ",CHI02,":",CMI02,":",CSI02))'

Where format included:
  * C(...) : print in calendar format 
  * CDI02 : calendar format for day as 2 digit integer with leading zero
  * CMoI02 : calendar format for Month as 2 digit integer with leading zero
  * CYI04 : calendar format for Year as 4 digit integer with leading zeroes
  * CHI02, CMI02, CSI02 : calendar format for hour, minute, second as 2 digit integers with leading zeroes

It is possible to print seconds as float:

  PRINT, t, FORMAT='(C(CHI02,":",CMI02,":",CSf05.2))'


General syntax of a Format specifier is 
  [<repeat>]<Format code>[+|-]<width>[.<precision>[E<exponent>]]
  
Where <repeat> <width> and <precision> are integer numbers. 
  * <repeat> says that the format code is repeatedly used. 
  * <width> is the width of the whole expression, <width>=0 means that the width is adapted to the number,  with <width>=03 leading zeros are set to fill a three digit field.
  * <precision> for floating point numbers gives the digits after the point, with integers the resulting number is filled to this width with leading Zeros.
  * E<exponent> gives the digits of the exponent for scientific notation
  * Format specifiers can be grouped with brackets which may have also a <repeat> field in front.

<Format code> is a usually a single letter:

| <Format code> | meaning | \\
|      A        | alpha numerical string |\\
|      I        | Integer: decimal representation |\\
|      B        | Integer: binary notation |\\
|      O        | Integer: octal notation |\\
|      Z        | Integer: hexadecimal notation |\\
|   F or D      | Floating point number |\\
|   e or E      | Floating point number in Scientific notation with 'e' or 'E' for the exponent |\\
|   g or G      | Floating point number in Scientific notation only if necessary |\\
|      C        | 'Calendar data' see [[idl:idl_jan#print formatted time|'julian day' below]] |


There are further formatting codes -> type '? format codes' at the idl prompt. Or look at [[https://www.nv5geospatialsoftware.com/docs/using_explicitly_formatt.html|using explicitly format]], [[https://www.nv5geospatialsoftware.com/docs/format_codes_fortran.html|format codes fortran]] or [[https://www.nv5geospatialsoftware.com/docs/format_codes_cprintf.html|format codes cprintf]].


==== mean, moments ====

function %%mean%% calculates the average of a complete array. 
  xm = mean( x , /nan )
Flag /nan makes idl to ignore nans.
If you also need variance and higher moments use instead moment which does the calculation faster:
  m = moment( x, /nan )
  x_mean = m[0]
  x_var  = m[1]
  s_skew = m[2]
  s_kurt = m[3]

With version 8 idl can restrict moment and mean to (only) one dimension:
  ; generate some data
  N = 1000
  x = randomu(rs,N) ; this gives a fltarr(N)
  y = (randomu(rs,N))^2; this gives a fltarr(N)
  xy = [ [x], [y] ] ; this gives a fltarr(N,2)
  m_x_y = mean(   xy, dimension = 1 ) ; this gives a fltarr(2) equal to [ mean(x), mean(y) ]
  m_xy  = mean(   xy, dimension = 2 ) ; this gives a fltarr(N) equal to m_xy[i] = mean([ x[i], y[i] ])
  mom = moment(   xy, dimension = 1 ) ; this gives a fltarr(2,4) equal to  [ [moment(x)], [moment(y)] ]
with variable dimension indicating the number of the dimension over which to average. The first dimension has number 1.

The same can be done with the max and min functions. 
  max_x_y = max( xy, dimension=1 ) ; this gives a fltarr(2) equal to [ max(x), max(y) ]
  max_xy  = max( xy, dimension=2 ) ; this gives a fltarr(N) equal to max_xy[i] = max([x[i], y[i]]) 
  ...

===== speed =====

IDL becomes slow if you use loops like for, while or repeat and conditional statments like if or case and even slower if you combine them.

This webpage provides some tips to speed your code: [[http://www.danidl.co.uk/idl.speed.tips.shtml]]

You can speed up your code if you use idl routines, matrix operations etc. Some examples below:


==== running average ====

We want a running average over 2*M points:

this is slow
  M = 11
  for i = 0, N-1 do begin
    ia = max([i-M,0])
    ib = min([i+M,N-1])
    xm[i] = mean( x[ia:ib] , /nan )
  endfor
  
use instead
  xm = smooth( x, 2*M , /edge_truncate , /nan )
The /nan key word let idl ignore NaN((not a number))

or in two dimensions if you want to average in x direction only.
  Mx = 2*M
  My = 1
  xm = smooth( x, [Mx,My] , /edge_truncate , /nan )

alternatively:
  Mx = 2*M
  My =  1
  kern = make_array( Mx, My, /float, value = 1.0/(Mx*My)  ) ; equal to kern=fltarr(Mx,My) & kern[*]=1.0/(Mx*My)
  xm = convol( x , kern , /edge_truncate , /nan )
Adapt Mx and My if you want more ...
This opens also possibilities to more sophisiticated smoothing techniques. Eg. with a gaussian smooother equivalent to the arithmetic mean over M elements:
  ; sigma derived from arithmetic average width M
  sigma_gauss = M/3.
  ; width should be at least 3sigma ... we take here 4 sigma
  N_gauss = floor(sigma_gauss*4)*2+1
  ; kernel is gaussian bell 
  kern = exp( - 0.5 * ( (findgen(N_gauss)-0.5*(N_gauss-1)) / sigma_gauss )^2 ) 
  ; nomrlaize it
  kern = kern/total(kern)
  ; convolute
  xm = convol( x , kern ,  /edge_truncate , /nan )

a running standard deviation is also easily done
  x_sdev = sqrt( smooth( (x-xm)^2 , 2*M , /edge_truncate , /nan ) )


==== variances, covariances ====

Instead of for-loops use the [[http://www.physics.nyu.edu/grierlab/idl_html_help/V2.html|variance]] or [[http://www.physics.nyu.edu/grierlab/idl_html_help/M39.html|moment]] routines of idl.

  var_x = variance(x , /nan )
  ; or 
  moms_x = moment( x , /nan  )
  mean_x = moms_x[0]
  var_x = moms_x[1]

For the **cross covariance** use

  xm = mean(x,/nan)
  ym = mean(x,/nan)
  covar_xy = mean( (x-xm)*(y-ym) , /nan )

If your data set is small you should multiply a factor N_data/(N_data-2) where %%N_data = total(finite(x-xm) AND finite(y-ym))%%

for the **time laged cross covariance** use:

  N_data = n_elements(x)
  xp = x - mean(x,/nan)
  yp = y - mean(y,/nan)
  N_lag = 10 ; must be smaller than N_data !
  covar_xy_pos = fltarr(N_lag) ; for positive lags = indgen(N_lag)
  covar_xy_neg = fltarr(N_lag) ; for negative lags =-indgen(N_lag)
  for i_lag = 0, N_lag-1 do begin
    covar_xy_pos[i_lag] = mean( xp[i_lag:N_data-1      ] * yp[  0  :N_data-1-i_lag] , /nan )
    covar_xy_neg[i_lag] = mean( xp[  0  :N_data-1-i_lag] * yp[i_lag:N_data-1      ] , /nan )
    endfor ; i_lag
  
And similar for the auto covariance by replacing y by x.


If you want to evaluate it for N_lag in the order of N_data it is much more time efficient to use the identity((see e.g. [[http://en.wikipedia.org/wiki/Cross-correlation#Properties|wikipedia]] about cross correlation)): 
  cov(x,y) = FFT( FFT(x) * CONJ(FFT(y)), /INVERSE )
but then there must be *not any NaN* in your dataset !\\
Be aware that:
  * covariances for negative timelags are at the end of the iFFT array.\\
  * this identity assumes periodicity of the time series i.e. x[i] = x[i+N], ... etc.



==== conditional if's ====

If you want to use data only if it is above a threshold this will be slow:
  for i = 0, N-1 do begin
    if x[i] lt threshold then begin
      x[i] = !values.f_nan
    endif ; x[i]
  endfor ; i

use instead:
  i_ts = where( x lt threshold , N_ts )
  if N_ts gt 0 then x[i_ts] = !values.f_nan
the if N_ts gt 0 is necessary because i_ts will be equal to -1 if no data is found which will result in an error prior to version IDL 8. But IDL 8 treats negative indices as relative indices - so you might not even notice that something wrong happened.

And similar if you do not want to juggle with NaNs create a new dataset:

  i_good = where( x gt threshold, N_good )
  if N_good gt 0 then x_good = x[i_good] $
                 else x_good = [] ; this gives !NULL (a pointer pointing nowhere ..) which may result in strange results


===== Programs and Routines =====

have the basic structure:

  PRO <name> [,<p1> [,<p2> [,<p3> [...]]]] [,<v1>=<v1> [,<v2>=<v2> [,<v3>=<v3> [...]]]]
  ...
  command block 1
  ...
  END
  [...
  command block 2
  ...
  END  ]

  * <name>   name of the program, it must be identical with the filename, except for the extension.\\ The code above must go into the file <name>.pro
  * <p1>,<p2>,... names of parameters
  * <v1>=<v1>,<v2>=... names of named parameters
  * command block 1  will be executed if the program is called
  * command block 2  is optional and will be executed if the program is compiled...

Example:
  PRO test ,x ,y ,dings=dings
  
  IF dings THEN PRINT,x,y
  END

must be saved as test.pro\\
and be compiled with 
  .RUN test
it can then be called with
  test,'aha',3.14,/dings
  test,'aha',3.14,dings=0
  test,'aha',3.14,dings=1
  test,'aha',/dings
  ...
some of these examples will produce errors...

To make parameters optional you should test whether they are present when the programm is called.
This can be done with:
  IF NOT ARG_PRESENT(dings) THEN dings = 0
ARG_PRESENT checks whether the variable dings was passed to the program\\
Nearly the same can be achieved with KEYWORD_SET which checks whether a variable **exists** and whether its **value is not zero**
  IF NOT KEYWORD_SET(dings) THEN dings = 0

===== standalone executables =====

We have only a limited number of licenses for IDL (60 but every user needs 6 at once thus 10effective, of these are 6 for idl8.2 and 4 for 6.4).
To reduce the time that IDL is blocked you may either build a standalone executable or use batch jobs:

==== standalone executables ====

In principle it is possible to generate standalone executables. They run in IDL virtual machine which executes the code
((so they are not really standalone)). Additional the virtual machine always asks for a mouse click - also if your application has no graphical user interface. The advantage is that it does not block a license.

Details how to generate and execute standalones are describe under "The idl virtual machine": \\
[[http://www.physics.nyu.edu/grierlab/idl_html_help/distributing6.html#wp1004661]]

And a more compact description at idl coyote: \\
[[https://www.idlcoyote.com/tips/idl_icon.html]]

Here is how to do it:
  * assume you have a main program **main.pro** which you want to run. In IDL execute the following command:
<code>
    .COMPILE main  ; compile main program
    RESOLVE_ALL    ; resolve all dependencies ...
    ; RESOLVE_ALL, CLASS=<object_names> ; ... if you use objects provide their names as a strarr in <object_names> 
    SAVE, /Routines, Filename='main.sav' ; generate the ''.sav'' file which may run in the virtual machine
</code>
  * in windows you may just click on **main.sav** to run it
  * in linux call from the command prompt:
<code>
    $ idl -vm=main.sav
</code>


==== batch jobs ====

If you want to run the idl from a bash script you may use an IDL batch job.

An idl program which name starts with a '@' character is interpreted as a batch job. It can be started from the comandline and occupies an idl license only as long it is running. You may edit the code between calls without blocking IDL.

Here is how to do it:
  * assume your idl stuff is done in **idl_program.pro** 
  * write a batch file '@call_idl_program'  (you can choose any name, but it must start with a '@' character)

  ; @call_idl_program
  ; compile and call idl_program.pro
  .COMPILE idl_program ; compile your stuff
  idl_program          ; run it
  EXIT                 ; leave idl

  * call it from command prompt with 

  $ idl @call_idl_program
  
  
You cannot pass paramters from the command line to your program. You either have to 
  * provide paramters in your  @call... file, e.g. generate it from your calling script
  * or call an idl program from within the @call... batch job
  * or provide parameters via an input file you read in idl_program.pro

a simple idl program to provide paramters could look like
  a = 1
  time = julday( 1,1,1970, 0,0,0 )
  x = 3.14*a
  end
save as set_my_params.pro and in @call_idl_program just include a line 
  ; @call_idl_program
  ; compile and call idl_program.pro
  .COMPILE idl_program ; compile your stuff
  set_my_params        ; set parameters
  idl_program          ; call your program
  EXIT                 ; leave idl

===== Functions =====

have in principle the same structure as programs, except that they return a value:

  FUNCTION <name> [,<p1> [,<p2> [,<p3> [...]]]] [,<v1>=<v1> [,<v2>=<v2> [,<v3>=<v3> [...]]]]
  ...
  command block
  ...
  RETURN, <variable>
  
  END

Example:
  FUNCTION sinx,x
  RETURN sin(x)/x
  END
===== Files  =====
==== text  ====

open (text)files for **input** into idl:

  s = '' ; this line predefines s as a string and the file is read textline by textline
  
  OPENR, file_unit, filename, /Get_LUN
  ...
    ; read one textline
    READF, file_unit, s
  ...
  CLOSE, file_unit
  FREE_LUN, file_unit

if the file is compressed with gzip it can be directly read by setting the flag %%/compress%% in %%OPENR%%

open (text)files for **output** from idl:

  OPENW, file_unit, filename, /Get_LUN
  ...
    PRINTF, file_unit, s
  ...
  CLOSE, file_unit
  FREE_LUN, file_unit


Read a text file containing a table with known format

  OPENR, file_unit, filename, /Get_LUN
  ...
    READF, file_unit, var_int, var_float, var_sci_float, format='(I3, F7.2, E12.2)'
  ...
  CLOSE, file_unit
  FREE_LUN, file_unit

Read a text file containing a table with columns separated by TAB's
<code idl>
  TAB = string(9b)
  OPENR, file_unit, filename, /Get_LUN
  ...
    READF, file_unit, s
    col = STRSPLIT( s , TAB , count=N_col , /extract )
  ...
  CLOSE, file_unit
  FREE_LUN, file_unit
  
</code>

Count **lines in a text file** can be done with function %%file_lines, filename%%

==== read_ascii ====

The above are simple alternatives to the very complex idl routine [[ http://www.harrisgeospatial.com/docs/READ_ASCII.html | READ_ASCII ]]:

  data = read_ascii( filename )

will read the data into a strucutre containing one array:

  help, data
  Structure <1485998>, 1 tags, length=<N_col>*<N_row>*<size_of_float>, data length=<N_col>*<N_row>*<size_of_float>, refs=1:
    FIELD1          FLOAT     Array[ <N_col>, <N_row> ]

it is possible to exclude comments, a header, or read only parts of the file by defining a %%TEMPLATE%%. 

This can be done with the function [[ http://www.harrisgeospatial.com/docs/ascii_template.html | ASCII_TEMPLATE]] which opens an interactive dialog.

  template = ASCII_TEMPLATE(  data_file_name )

In principle this function can detect several features of the file (like column separator, data type per column, ...)
But you have to click several buttons.

The function returns a structure with the following fields:

  Structure <179f0e8>, 10 tags, length=120, data length=117, refs=1:
   VERSION         FLOAT           1.00000
   DATASTART       LONG                 0
   DELIMITER       BYTE        <delimiter ascii code>
   MISSINGVALUE    FLOAT               NaN
   COMMENTSYMBOL   STRING    '<comment symbol>'
   FIELDCOUNT      LONG         <N_col>
   FIELDTYPES      LONG      Array[<N_col>]
   FIELDNAMES      STRING    Array[<N_col>]
   FIELDLOCATIONS  LONG      Array[<N_col>]
   FIELDGROUPS     LONG      Array[<N_col>]

Accordingly it should be possible to construct it in the code to avoid calling the dialog:

<code>

; delimiter between columns
; delimiter = ',' ; comma for csv files
delimiter = string(9b) ; tab for tsv files

; read header
; the line where to find the column header , we assume here the very first line
i_header = 0
header = '' ; predefine as a string => file can be read textline by textline
OPENR, file_unit, filename, /Get_LUN
; read first textlines until we have the column header
for i = 0, i_header do READF, file_unit, header
;   one may alternatively search for a certain key in the text line with
;   key = '# time'+delimeter
;   N_header = 0
;   WHILE (strpos(header,key) lt 0) and (~ eof(file_unit)) DO BEGIN
;     READF, file_unit, header
;     N_header += 1
;     ENDWHILE ; key in header
;   i_header = N_header-1
; check whether there remains anything in the file
if eof(file_unit) then message, 'no data in file !'
; close file
CLOSE, file_unit
FREE_LUN, file_unit

; split header, and count columns separated by delimiter.
header = strsplit( header, delimiter, /extract, /preserve_null, count = N_col )

; define type of variables in different columns
; type: 0=skip field, 1=byte, 2=int16, 3=long32, 4=single, 5=double, 7=string, ... see idl help for size
; i assume here that the first 3 columns are strings, i.e. cannot be converted into numeric values and the rest can be converted to single precision floating point
N_col_str = 3
fieldtypes = [ replicate(7,N_col_str), replicate(4,N_col-N_col_str)  ]

template = { $
  VERSION         : 1.0, $ ; must be 1.0 ! IDL 6.0 does not know anthing else ...
  DATASTART       : i_header+1 , $  ; index to line where to start reading. First line=0. We assume here that data starts directly after the header
  DELIMITER       : delimiter, $ ; one character - must be defined above e.g. as TAB=string(9b) or comma=',' ...
  MISSINGVALUE    : !values.f_nan, $ ; idl inserts this if it cannot read
  COMMENTSYMBOL   : '#', $ ; one character ...
  FIELDCOUNT      : N_col, $ ; number of columns, number of elements of following fields must fit
  FIELDTYPES      : fieldtypes, $ ; type: 1=byte, 2=int, ... 4=single, 5=double, 7=string, see idl help for size
  FIELDNAMES      : header, $ ; these names will be used for the variables, i.e. entries in header must follow the rules for variable names
  FIELDLOCATIONS  : replicate(  0, N_col ), $ ; indices to column begin, if no delimiter is given (?), or if columns are of constant width
  FIELDGROUPS     : indgen(N_col) $ ; every field an independent group, i.e. one element in the struct data we are going to read
  ; FIELDGROUPS     : replicate(0,N_col) $ ; every field the same group - only possible if the same type ! => data = { header[0] : array(N_col:N_lines) }
  ; FIELDGROUPS     : [ replicate(0,N_col_str) , replicate(1,N_col-N_col_str) ] ,  $ ; one group for the strings, one for the floats => data = { header[0] : strarr(N_col_str,N_lines), header[N_col_str] : fltarr(N_col-N_col_str,N_lines) }
  }

; read the data
data = read_ascii( filename, template=template )

; show info about the structure of your data
help, data, /struct

...

</code>

  * FIELDCOUNT and number of elements of FIELDTYPES, -NAMES, -LOCATIONS and -GROUPS must be the same !
  * if you group your data 
    * you will get for every group an %%M x N%% array with M the number of group members and N the number of records (or lines)
    * you nevertheless need to provide FIELDTYPES, -NAMES, -LOCATIONS and -GROUPS with FIELDCOUNT elements !
    * the name of every group is then the name of the first FIELDTYPENAME going into it


==== binary  ====

with READU **binary data** can be read: 
  ...
  ; define x as long variable
  x = 0l
  ; read it from file
  READU, file_unit, x
  ...

with WRITEU binary data can be written:
  ...
  WRITEU, file_unit, x
  ...


==== file size, position etc ==== 

check for existence
  file_info_rec = FILE_INFO( filename )
  if file_info_rec.exists then print,'exists' $
                          else print, 'does not exist'
                          
or in one line
  if (FILE_INFO( filename )).exists then print,'exists' else print, 'does not exist'

get file size 
  file_info_rec = FILE_INFO( filename )
  print, 'file size=', file_info_rec.size
  

set file pointer to a certain position
  OPENR, file, filename, /Get_LUN
  ...
  POINT_LUN, file , file_position

get position of file pointer (=> negative file unit number!) :
  POINT_LUN, -file , file_position

count lines in a text file
  N_lines = FILE_LINES( file_name ) 


===== NETCDF =====

We prefer to save data as **NETCDF files** (=**N**etwork **C**ommon **D**ata **F**ormat) \\ 
NetCDF is a platform independent binary format which includes meta information about its content. \\
It can be easily read in many programming platforms like matlab, python, C, ... \\
As a binary format reading is very fast and access to a subset of the content is easy.\\

  * NETCDF was developed by [[http://www.unidata.ucar.edu/software/netcdf/|Unidata]]@[[http://www.ucar.edu|UCAR]] see also [[http://de.wikipedia.org/wiki/NetCDF|en.wikipedia.org]] ...
  * a description can be found at [[https://docs.unidata.ucar.edu/netcdf-c/current/index.html]]
  * beside the values they should contain **meta-information** called **attributes** \\ this can be units of the values, instrument location, model run number  etc.\\ unidata provides [[http://www.unidata.ucar.edu/software/netcdf/docs/netcdf/Attribute-Conventions.html#Attribute-Conventions|conventions about attributes]]
  * the data is organised in arrays i.e. every data record in a netcdf file has at least one **dimension**
  * IDL netcdf overview: [[https://www.l3harrisgeospatial.com/docs/NCDF_Overview.html]]
  * NetCDF error codes [[https://docs.unidata.ucar.edu/netcdf/documentation/4.8.0/nc-error-codes.html]]


I wrote functions **%%nc_read.pro%%** and **%%read_nc_group.pro%%** to read netcdf files into a structure comparable to the Python function [[https://docs.xarray.dev/en/stable/generated/xarray.open_dataset.html#xarray.open_dataset|xarray.open_dataset]]. \\
You can find them in %%/home/jschween/exp/ideen/idl/netcdf%%.

For details of IDLs netcdf functions see the IDL help pages for a [[https://www.nv5geospatialsoftware.com/docs/NCDF_Overview.html|NCDF Overview]].\\
Below example codes rely on the basic NETCDF functions IDL provided since long. Since version 8.4.1 IDL provides a //simplified interface// which also provides compact access to netcdf data.

==== Create a NETCDF file ====

  ; create netcdf output file
  ; use flag /clobber => overwrite exisiting file (the verb to clobber means to hit or punish someone, what strange... )
  nc_file = NCDF_CREATE( nc_filename , /clobber )
  
  ; set global attributes
  NCDF_ATTPUT, nc_file, /Global, 'Creator' , 'Jan Schween (jschween@uni-koeln.de)'
  NCDF_ATTPUT, nc_file, /Global, 'creation_date' , STRING( SYSTIME(/utc,/jul), f='(C(CDI02,".",CMoI02,".",CYI04,"/",CHI02,":",CMI02))' )
  
  ; define time as dimension
  ; the corresponding variable shall have the same name as it is intended to be a coordinate variable
  N_time_id = NCDF_DIMDEF( nc_file, 'time', N_time )
  
  ; we use NaN for missing values
  novalue = !values.f_nan
  
  ; define time variable with dimension time, 
  ;   as time is a coordinate variable its name should be the same as the dimension
  ;   additionally as coordinate variable it should be in ascending order and have no gaps (missing values, filled values)
  ;   we use here the IDL - time format julian day which allows the determination of calendrial dates or can be inversely determined from calendar dates.
  time_id = NCDF_VARDEF( nc_file, 'time', N_time_id, /DOUBLE )
  ; set attribute for the units of time
  NCDF_ATTPUT, nc_file, time_id, 'units' , 'days since -4712-1-1T12:0:0' ; format according to cfconventions, see https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html , this should work with   "ncdump -iv time <nc_file>"   ... not allways but sometimes(?) ...  
  NCDF_ATTPUT, nc_file, time_id, 'calendar', 'gregorian' ; this is the standard calendar - in principle you can omit this.
  NCDF_ATTPUT, nc_file, time_id, 'long_name', 'time as julian day'
  NCDF_ATTPUT, nc_file, time_id, 'missing_value', double(novalue)
  NCDF_ATTPUT, nc_file, time_id, '_fillvalue', double(novalue) ; allthough unidata gives "_FillValue" with capital letters, you may get an error when using them... (see [[http://www.rhinocerus.net/forum/lang-idl-pvwave/575253-ncdf_attput-_fillvalue-problem-string-arrays.html|here]])
  NCDF_ATTPUT, nc_file, time_id, 'comment', 'Julian day are the fractional days since 1.1.4713BC 12GMT. It can be converted to unix epoch (t_unix:seconds since 1.1.1970) with t_unix=(time-2440587.5)*86400.'
  
  ; alternativley time as unix epoch - wide spread time format usually of type integer, but then with overflow problem...
  tunx_id = NCDF_VARDEF( nc_file, 'tunix', N_time_id, /LONG ) 
  NCDF_ATTPUT, nc_file, tunx_id , 'units' , 'seconds since 1970/1/1 0:0:0'  ; format according to cfconventions, , evtly. '-' as separator in date? see https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html, no time zone info means automatically UTC
  
  ; create id for a variable dependent on time, 
  ; possible types are: /BYTE , /CHAR, /DOUBLE , /FLOAT , /LONG , /SHORT, /STRING, /UBYTE, /UINT64, /ULONG, /USHORT
  ; where /SHORT is a two byte integer, /LONG is a 4byte integer, and /UINT64 is an unsigned 8byte integer,
  ; and you may already guess: /U... means that the variable is an unsigned integer
  var_id = NCDF_VARDEF( nc_file, 'some variable name', N_time_id, /FLOAT )
  ; multidimensional arrays like fltarr(N_time,N_space) are defined via
  var_2d_id = NCDF_VARDEF( nc_file, '2D variable', [ N_time_id, N_space_id ], /FLOAT )
  ; single vars are defined by omitting the dimension
  var_one_value_id = NCDF_VARDEF( nc_file, 'one_value', /FLOAT )
  
  ; set attributes
  ; these here are standard attributes (including the underscore) as defined in the unidata conventions ...
  NCDF_ATTPUT, nc_file, var_id, 'long_name', 'long name'
  NCDF_ATTPUT, nc_file, var_id, 'units', 'arb.units'
  NCDF_ATTPUT, nc_file, var_id, 'missing_value', novalue
  NCDF_ATTPUT, nc_file, var_id, '_fillvalue', novalue ; allthough unidata gives "_FillValue" with capital values, you may get an error when using them... (see [[http://www.rhinocerus.net/forum/lang-idl-pvwave/575253-ncdf_attput-_fillvalue-problem-string-arrays.html|here]])
  
  
  
  ; further variables and attributes ...
  ; ...
  
  ; STRINGS are stored as array of (unsigned) bytes called CHAR
  ; we therefore need the length of the string
  ; and the netcdf vaiable has to dimensions:
  ; lets assume we have str_var = strarr(N_str_var) 
  max_str_len = max( strlen(str_var) )
  ; define dimension
  str_len_id = NCDF_DIMDEF( nc_file, 'str_len', max_str_len )
  str_var_id = NCDF_VARDEF( nc_file, 'str_var_name', [ str_len_id, N_str_var_id ], /CHAR )
  
  ; close define mode
  NCDF_CONTROL, nc_file, /endef
  
  ; put variable into file ...
  NCDF_VARPUT, nc_file, time_id, time
  NCDF_VARPUT, nc_file, tunx_id, (time-2440587.5)*86400.
  
  ; STRING variable is as simple as the others ...
  NCDF_VARPUT, nc_file, str_var_id, str_var
    
  ; close netcdf file
  NCDF_CLOSE, nc_file
  

==== NETCDF time.units ====

For time.units I follow here [[https://cfconventions.org/Data/cf-conventions/cf-conventions-1.7/build/ch04s04.html|cfconventions]].\\
But they are not very clear - they just give examples which differ from each ohter.\\
Accordingly time.units could have the form:
  * '-' as separator in date
    * ''<base_unit> since <[-]YYYY-MM-DD hh:mm:ss>''
  * clock time with no leading zeros
    * ''<base_unit> since <[-]YYYY-MM-DD [h]h:[m]m:[s]s>''
  * '/' as separator in date
    * ''<base_unit> since <[-]YYYY/MM/DD hh:mm:ss>''

Commandline tool ''ncdump'' cannot deal in every case with the clock-time in time.units.\\
what seem to work with ncdump is:
  * you need dashes '-' as separators in date 
  * blank ' ' or 'T' as separator between date and time
  * Leading zeros can be omitted
  * that means
    * ''[-][Y[Y[Y]]]Y-[M]M-[D]D{ |T}[h]h:[m]m:ss''



==== Dimensions ==== 

In the above example 'time' is the name of the **dimension** and of a **variable**. That is intended to make clear the close relation between %%N_time%% and %%time = dblarr(N_time)%% and part of the [[http://cf-pcmdi.llnl.gov/|CF netcdf dataformat]]((CF=Climate and forecast)). But it is also possible to use different names, in this case ncview allows you to plot an array against other variables...

In netcdf 64 it is possible to have one 'UNLIMITED' dimension. I.E. this dimension does not need to be specified at creation time and data can be appended to the file. In idl this is done by using the keyword '/unlimited' (see also the example 'appending data ...' below):
  id_dim_time = NCDF_DIMDEF( nc_file, 'time', /unlimited )



==== Attributes ==== 

There is a list of 
[[http://www.unidata.ucar.edu/software/netcdf/docs/netcdf/Attribute-Conventions.html#Attribute-Conventions|netcdf attribute conventions]]
at the unidata pages


==== possible datatypes in netcdf: ==== 

IDL help is not 100% clear about the meaning of the datatypes in netcdf :

  * %%/BYTE%% 8bit signed/unsigned(?) integer ... ncview assumes signed i.e. numbers between -128 and 127
  * %%/CHAR%% this is also byte, but //assumed to be characters// i.e. a string can be stored as an array of char but you have to convert it in idl to string. Might be also interpreted as unsigned 8bit integers i.e. numbers between 0 and 255
  * %%/SHORT%% 16bit integer (2 byte) probably signed, i.e. numbers between -32768 and +32767
  * %%/LONG%% //longword// integers i.e. probably 32 bit signed integers i.e. numbers from − 2,147,483,648 to 2,147,483,647
  * %%/FLOAT%% floating point numbers - probably [[http://en.wikipedia.org/wiki/Single_precision_floating-point_format | IEEE single precision]] (32bit)
  * %%/DOUBLE%% [[http://en.wikipedia.org/wiki/Double_precision_floating-point_format | IEEE double precision]]  (64bit)


==== Read a NETCDF file:====

  ; open netcdf file
  nc_file = NCDF_OPEN( filename , /NOWRITE )
  
  ; get global attribute
  NCDF_ATTGET, nc_file, /GLOBAL, att_name, glob_att
  
  ; get id of dimension
  dim_id = ncdf_dimid( nc_file , 'dim_name' )
  
  ; get id of variable via its name
  var_id = NCDF_VARID( nc_file, 'some variable name' )
  
  ; get attribute for variable
  NCDF_ATTGET, nc_file, var_id, 'units', var_unit_str
  
  ; get dimension
  NCDF_diminq, nc_file , dim_id, name_dim , N_dim 
  
  ; elements of a metadata group
  ; see on Harris
  ;   http://www.harrisgeospatial.com/docs/ncdf_ncidinq.html
  ;   http://www.harrisgeospatial.com/docs/ncdf_groupname.html#netCDF_2618656010_1006782
  metadata_grp_id = NCDF_NCIDINQ( nc_file , Metadata_group_name )
  mata_data_id = NCDF_VARID( metadata_grp_id, meta_data_name )
  
  ; get data
  NCDF_VARGET, nc_file, var_id, variable
  
  ; get string data
  str_var_id = NCDF_VARID( nc_file, 'str_var_name' )
  NCDF_VARGET, nc_file, str_var_id, str_var
  str_var = string( str_var )
  
  ; close netcdf file
  NCDF_CLOSE, nc_file


==== appending data to an exisitng netcdf file: ====

Appending data is possible if the dimension is defined as //unlimited//. 
Defining the variables is as usual with one exception:

  ; define unlimited dimension time 
  id_dim_time = NCDF_DIMDEF( nc_file, 'time', /unlimited ) ; only one inlimited dimension allowed (?)
  id_dim_2D = NCDF_DIMDEF( nc_file, 'N_2D', N_2d )
  ; connect variables to it
  id_time = NCDF_VARDEF( nc_file, 'time', id_dim_time , /DOUBLE)
  id_data = NCDF_VARDEF( nc_file, 'data', id_dim_time , /FLOAT )
  id_data_2D = NCDF_VARDEF( nc_file, 'data', [ id_dim_2D, id_dim_time ], /FLOAT ) ; for 2D data the unlimited dimension must be at the end

If the file remains open just append the data at the desired position

  ; append data assuming that variables //time// and //data// contain only new data
  NCDF_VARPUT, nc_file, id_time, time, offset=[ N_time ]
  NCDF_VARPUT, nc_file, id_data, data, offset=[ N_time ]
  NCDF_VARPUT, nc_file, id_data_2D, data, offset=[ 0, N_time ]

If the file is reopened the number of elements already in the file must be determined:

  nc_file = NCDF_OPEN( nc_file_name , /write )
  ; get id for variables time and data
  id_time = NCDF_VARID( nc_file , 'time' )
  id_data = NCDF_VARID( nc_file , 'data' )
  
  ; get number of elements in time in two steps
  info_time = NCDF_VARINQ( nc_file , id_time )
  NCDF_DIMINQ, nc_file, info_time.dim[0], name, N_time
  
  ; alternatively
  id_time = NCDF_DIMID( nc_file , 'time' )
  NCDF_DIMINQ, nc_file, id_time, name, N_time
  
  ; append data assuming that variables //time// and //data// contain only new data
  NCDF_VARPUT, nc_file, id_time, time, offset=[ N_time ]
  NCDF_VARPUT, nc_file, id_data, data, offset=[ N_time ]


==== Compression ====

With NETCDF 4 it is possible to compress data (see [[https://www.unidata.ucar.edu/blogs/developer/entry/netcdf_compression|unidata.ucar.edu]]). 
IDL can read NETCDF4 file as of version 8.0 (see IDL ref pages for [[https://www.nv5geospatialsoftware.com/docs/ncdf_create.html|NCDF_CREATE]]).
You can check the NETCDF version implemented in IDL with %%HELP, 'ncdf', /DLM%%.

With our IDL9 implementation this does not work anymore: you get an netcdf error -31 after NCDF_CREATE !?!?\\
Although it say it has version 4 implemented:

  IDL> help, 'ncdf' , /DLM
  ** NCDF - IDL Network Common Data Format (NetCDF) support DLM (loaded)
    Version: netcdf-c 4.9.2, Build Date: NOV  8 2023, Source: NV5 Geospatial Solutions, Inc.
    Path: /opt/nv5/idl90/bin/bin.linux.x86_64/idl_netcdf.so


To use NETCDF compression you need to set the NETCDF4_FORMAT keyword when generating the file,\\
and the GZIP and SHUFFLE keywords when defining a variable, eventually also CHUNK_DIMENSIONS :

  ; create netcdf output file
  nc_file = NCDF_CREATE( filename , /clobber, /NETCDF4_FORMAT )
  ...
  g_zip_level = 5 ; a value between 0..9 for low...high compression but also fast...slow access
  ; chunk dimensions: can help the compression but can also be omitted
  cd_time  = N_time / 4 ; i guess that is the way to use it
  cd_range = N_range / 4 ; i guess that is the way to use it
  
  ; get id for a variable dependent on time, 
  ; possible types are: /BYTE , /CHAR, /DOUBLE , /FLOAT , /LONG , /SHORT
  data_id = NCDF_VARDEF( nc_file, 'data', N_time_id, /FLOAT , GZIP=g_zip_level, CHUNK_DIMENSIONS=cd_1, /SHUFFLE )
  ; multidimensional arrays are defined via
  data_2d_id = NCDF_VARDEF( nc_file, 'data_2d', [ N_time_id, N_range_id ], /FLOAT , GZIP=g_zip_level , CHUNK_DIMENSIONS=[cd_time, cd_range ], /SHUFFLE )
  ...

  * GZIP determines the compression rate. Allowed values are 1..9 for lowest to highest compression or fastest to slowest compression (see [[https://www.nv5geospatialsoftware.com/docs/ncdf_vardef.html#netCDF_2618656010_1007336|help for NC_VARDEF#GZIP]]).
  * CHUNK_DIMENSIONS determines how large blocks are compressed, it is by default the whole size of the variable except for variable of undefined length (see [[https://www.nv5geospatialsoftware.com/docs/ncdf_vardef.html#netCDF_2618656010_1007314|help for NC_VARDEF#CHUNK_DIMENSIONS]]). This means you can omit it. I dont know whether this paramter can improve compression.
  * SHUFFLE tries to sort bytes before compression (see [[https://www.nv5geospatialsoftware.com/docs/ncdf_vardef.html#netCDF_2618656010_1007349|help for NC_VARDEF#SHUFFLE]]).
  * with just GZIP=9 i achieve compression ratios around only Size_comp/Size_org = 85% for some smooth data, when appliying /SHUFFLE this improves to ~40%


==== Netcdf-groups ====

Data in Netcdf files can be organized in groups.
  * to get a list of the ID's  all groups in a netcdf file use [[https://www.nv5geospatialsoftware.com/docs/ncdf_groupsinq.html | NCDF_GROUPSINQ ]]  
  * to get the names of the groups use [[https://www.nv5geospatialsoftware.com/docs/ncdf_groupname.html#netCDF_2618656010_1006782 | NCDF_groupname ]]
  * to get the ID of a group from its name use [[https://www.nv5geospatialsoftware.com/docs/NCDF_NCIDINQ.html | NCDF_NCIDINQ ]]

Sample code:
  ; open file
  file_Id = NCDF_OPEN(fileName)
  ; get list of goup-id's
  group_id_list = NCDF_GROUPSINQ( file_id )
  for i = 0, n_elements( group_id_list )-1 do begin
    group_id_i = group_id_list[i]
    group_name_i = NCDF_GROUPNAME( group_id_i )
    print,  i, '->', group_id_i, '->', group_name_i
    ; in principle a group can again contain groups ...
    ; note that we now use the groups_id to reference the parent group where we used before the file id.
    sub_group_id_list = NCDF_GROUPSINQ( group_id_i )
    ... 
  endfor ; i


==== Rules for variable names in netcdf files: ====

One idea of netcdf is to read a file directly and automatically into variables with names equal to the names in the file.
That means that netcdf-variable names have to follow certain rules and not every character is allowed in a variable name in an netcdf file.

In general: a valid variable name in a programming language like //fortran// is an allowed variable name for a netcdf file. There is an [incomplete] description of the rules at the [[https://www.unidata.ucar.edu/software/netcdf-java/CDM/Identifiers.html| Unidata pages]]:

A valid identifier can be:

  ID = ([a-zA-Z0-9_]|{UTF8})([^\x00-\x1F\x7F/]|{UTF8})*

which says:
  * the first char must be alphanumeric, an underscore, or a UTF8 multibyte encoded char.
  * subsequent characters may be anything except control characters [\x00-\x1F] , forward slash [/], and del [\x7F]. So the excluded ascii decimal values are 0-31, 47, and 127.
  * A UTF-8 multibyte encoding always has the high order bit set ...

to be complete: the following are not allowed as well:
  * '(' and ')' 
  * '/' 
  * '%' 
  * white space
  * ':' (colon)
  * '"' (double quotes)
  * '*' '+' ... (operators)
  * etc. ...

The reason for this strict variable naming is that the variable names from inside the file should also be used as variable names in a program. 
Reading would be thus very simple: just read the file and you have the variables with the same name as in the file. 
Unfortunately this is not impelemented in idl...



===== Plot =====

  PLOT, <x-array>, <y-array>

generates a rectangle (box) with automatically scaled abscissa (bottom) and ordinate (left side) and plots  <y-array> versus <x-array>. 
There are many [[https://www.nv5geospatialsoftware.com/docs/PLOT_Procedure.html|optional parameters]] and many [[https://www.nv5geospatialsoftware.com/docs/graphkeyw.html#graphkeyw_3288778166_331699|addtional keywords]] to modify the plot.

The scaling is by default something like 10% wider than the range of the data but this can be modified by the [[https://www.nv5geospatialsoftware.com/docs/graphkeyw.html#graphkeyw_3288778166_332116|style]] keyword. You set essentially single bits of it and may combine (at least some of) them:

  [x|y|z]style
  1 -> Force exact axis range.
  2 -> Extend axis range.
  4 -> Suppress entire axis
  8 -> Suppress box style axis (i.e., draw axis on only one side of plot)
  16 -> Inhibit setting the Y axis minimum value to 0 (Y axis only) 

suppress entire x-Axis (top and bottom of the box) with style option
  xstyle=4
suppress top axis and right y-axis
  xstyle=8, ystyle=8
  


Suppress annotation of axis-ticks
  xtickname=replicate(' ',30)
Die Achseenbeschriftung besteht aus dreisssig Leerzeichen, IDL beschriftet die Achse nicht ...\\
... oder multiplot.pro von mir verwenden.

The rectangle around a plot is in idl called AXIS - there is a command AXIS to just define this rectangle without plotting it. But there is no way to return to a formerly defined one.

Contour plots (isolines and/or color shading ) can be done with the command [[http://www.physics.nyu.edu/grierlab/idl_html_help/C40.html#wp908085|countour]]:

  CONTOUR, z, x, y

I takes in principle the same options and keywords as the plot command plus some other.

Polar contours can be done with:

  POLAR_CONTOUR, field, az_list, R_list


===== ticks =====

you can define the length of intervals between ticks and the number of sub intervals with {x|y}tickinterval and {x|y}minor where minor gives the number of subintervals - not ticks !
If you want to have ticks every 2.5 units and minor ticks every  0.1 then:

  plot, ... , xtickinterval = 2.5, xminor=5

or 

  !X.xtickinterval = 2.5
  !X.minor = 5

if you do not want minor ticks set 

  xminor = 1
  


===== Annotation =====

you may want to annotate an axis with your own words.
You then have to set {x|y}ticks with the number of intervals (not ticks !) {x|y}tickname with an array of strings for the annotations and {x|y}tickv for the locations where to put the tics. Note that {x|y}ticks is the number of intervals i.e. n_elements({x|y}tickname)-1. If you omit xticks IDL assumes 5 intervals or 6 ticks !.

Here is an example with a logarithmic time axis:

  hours_per_day = 24d
  days_per_year = 365d
  
  ; basic unit of our x-axis is years - thus it is:
  hour = 1./(hours_per_day*days_per_year)
  day =  24*hour
  month = 30*day  ; ... well roughly - average length of a month is 30.416.. and 30.5days for a normal and a leap year, respectively, or ~30.44days on average
  
  x_tikvals = [ hour,   day,  month,   1.,   10. ,  100. ,  1000.,  10e3,   100e3,   1e6   ]
  x_tiknams = [ '1h',  '1d',  '1Mo',  '1y', '10y', '100y',  '1ky', '10ky', '100ky', '1My'  ]
  N_ticks = n_elements(x_tiknams)
  
  plot, [1], [1], $
    xstyle = 1, $ 
    xrange = [ 0.5 * hour, 2*1e6 ], $
    /xlog, $
    xticks = N_ticks-1, $  ; holy IDL - why do you need things like this ?
    xtickname = x_tiknams, $
    xtickv    = x_tikvals, $
    ystyle  = 1, $
    yrange = [ 0.5, 4.5 ], $
    yminor = 1, $
    /nodata

Similar for an x-axis based on julian day:

  hours_per_day = 24d
  
  ; basic unit of our x-axis is (julian)days - thus it is:
  hour = 1./hours_per_day
  day =  24*hour
  month = 30*day  ; ... roughly: average length of a month is 30.416.. and 30.5days for a normal and a leap year, respectively, or ~30.44days on average
  year = 365.*day ; ... roughly ...
  minute = hour/60.
  second = minute/60.
  milisec = second/1000.
  
  x_tikvals = [ milisec, second, minute,  hour ,  day ,  month,    year   ]
  x_tiknams = [ 'msec',   'sec',  'min', 'hour', 'day', 'month',  'year'  ]
  N_ticks = n_elements(x_tiknams)
  
  plot, [1], [1], $
    xstyle = 1, $ 
    xrange = [ 0.5*x_tikvals[0], 2*x_tikvals[N_ticks-1] ], $
    /xlog, $
    xticks = N_ticks-1, $  ; holy IDL - why do you need things like this ?
    xtickname = x_tiknams, $
    xtickv    = x_tikvals, $
    ystyle  = 1, $
    yrange = [ 0.5, 4.5 ], $
    yminor = 1, $
    /nodata

==== Annotation of a time axis ====

IDL provides in procedure PLOT the parameter XTICKUNITS which allows automatic annotation of a time axis. If t is time in julianday then

  plot, t, x, XTICKUNITS = 'Months'
  
generates a plot where every first of a month a tick is set. Annotations are the three first letters of the name of month (format = 'CMOA' )

Other possible values for XTICKUNITS are:

  * 'Days'  
  * 'Hours' 
  * 'Minutes'
  * 'Seconds'

A tick is set in fixed intervals optimal for the length of the axis (N days, or N hours etc.). 
Annotation is the day of the month or hour of day etc. but can be formatted with XTICKFORMAT.
For example:

  plot, t, x, XTICKUNITS = 'Days', XTICKFORMAT = '(C(CDI02,",",CMOI02))'

Will give annotations of the form DD.MO instead of the month name.

Another example:

  plot, t, x, XTICKUNITS = 'Hours', XTICKFORMAT = '(C(CHI02,":",CMI02))'

Will give annotations of the form HH:MM instead of just hour of the day as one or two digit number.


Another way is to explicitly define format and intervals, e.g. :

  plot , t, y , XTICKFORMAT = '(C(CHI02,".",CMI02,"."))', XTICKINTERVAL = 3/24., XMINOR = 3

annotates HH.MM every 3 hours (=3/24. of a day) and generates two minor ticks for 3 subintervals i.e. minor ticks at full hours.


You may want to have annotations in two lines: use code !C for a newline

  plot, t, y, XTICKFORMAT = '(C(CDI02,".",CMoI02,"!C",CYI04))'



==== xstyle: axis layout ====
  {x|y|z}style = 
   1 => force exact range
   2 => extended range
   4 => supress entire axis (dont draw line and tic marks ...)
   8 => supress box style - axis only one one side of plot
  16 => supress y=0 value for y axes - you can also use the /YNOZERO flag in the plot command...

==== linestyles ====

  linestyle = 
   0 => Solid
   1 => Dotted 
   2 => Dashed 
   3 => Dash Dot 
   4 => Dash Dot Dot 
   5 => Long Dashes 

==== symboltypes ====


  PSYM =
   1 => Plus sign (+) 
   2 => Asterisk (*) 
   3 => one pixel dot (.) 
   4 => Diamond 
   5 => Triangle up (if symsize>0)  or down if symsize<0
   6 => Square 
   7 => X 
   8 => User-defined. See USERSYM procedure. 
   9 => Undefined 
  10 => Histogram mode. Horizontal and vertical lines connect the plotted points, 
        as opposed to the normal method of connecting points with straight lines. 
        Vertical lines are located in the middle between subsequent x.
        See Histogram Mode for an example. 
        
Size of the symbol is controlled via the [[http://www.physics.nyu.edu/grierlab/idl_html_help/graphkeyw2.html#wp316061|symsize]] keyword. A value of one is about one charactersize.

IDL does not know filled symbols, but one can use alternatively the userdefined symbol PSYM=8 and the USERSYM procedure. 
To generate a filled triangle pointing to the top use:
<code idl>
  N_sym = 3
  w0_sym = 0.
  w_sym = findgen(N_sym+1)/N_sym*2*!Pi
  usersym, sin(w_sym+w0_sym) , cos(w_sym+w0_sym), /fill
  plot, x , y, psym=8
</code>
To make different symbols set:

|  symbol            |  N_sym  |  w0_sym  |
|  triangle 'up'     |  3      |    0     |
|  triangle 'right'  |  3      |   !pi/2  |
|  triangle 'down'   |  3      |   !pi    |
|  triangle 'left'   |  3      |  -!pi/2  |
|  diamond           |  4      |    0     |
|  square            |  4      |   !pi/4  |
|  circle            |  15     |    0     |

More complex symbols like stars might be generated by using something like

  usersym, sin(w_sym*N_arms)*sin(w_sym+w0_sym) , sin(w_sym*N_arms)*cos(w_sym+w0_sym), /fill
  
Or you use programs from the internet like [[http://idlastro.gsfc.nasa.gov/ftp/pro/plot/plotsym.pro|plotsym.pro]] from [[http://idlastro.gsfc.nasa.gov/]]


==== multiple data rows ====

To plot multiple data rows use the %%OPLOT%% command:
<code idl>
  ; plot y versus x
  PLOT, x, y
  ; plot y2 versus x on same axis
  OPLOT, x, y2
</code>
==== multiple axes ====

If you want to use several different axes use the %%AXIS%% command. To plot some data %%y%% and %%y2%% on the same x-axis but on the different y-axes on the (standard) left and the right side use the following code:
<code idl>
  ; define test data
  N = 1000
  x = findgen(N)/N
  y = sin(2*!pi*x)
  y2= 100.*exp(-((x-0.5)/0.1)^2)


  ; plot y data on (standard) left y-axis ...
  ; ... leave additional space of 10 charwidths for x-axis-labels on both sides
  PLOT, x, y1, ytitle='y1', XMARGIN=[10,10]
  ; define second y-axis on the right side , and save it for use
  AXIS, yaxis=1, yrange=[min(y2),max(y2)], ytitle='y2', /save
  ; plot data y2 onto this axis
  OPLOT, x, y2
</code>
The keyword %%yaxis=0%% refers to a y-axis on the left side of the plot with tics poiting to the right side (i.e. inward). The keyword %%yaxis=1%% refers to an y-axis on the right side of the plot with tics poiting to the left side (i.e. inward). Nevertheless the position of the axis can be set by giving the x-position of the axis. Here is an example for a third axis:

<code idl>
  ; additional data...
  y3 = atan(10.*(x-0.3))/(!pi/4.)
  ; position for third y-axis - left to the first axis
  x_icpt = -0.13
  ; define third y-axis placed at x=x_icpt, and save it for use
  AXIS, x_icpt, yaxis=0, yrange=[min(y3),max(y3)], ytitle='y3', /save
  ; plot data y3 onto this axis
  OPLOT, x, y3
</code>

Of course you can define additonal x-Axes you then have to use the parameter %%XAXIS=1%% (Xaxis on top) or %%XAXIS=0%% (x-axis at the bottom).



Sometimes the title/label for the second axis is not readable (especially if you are creating .ps-plots). An easy solution is to type: !p.position=[0.1,0.1,0.9,0.9] in front of the plot-commands.  (Stefan) ... More elaborated is to use the %%XMARGIN%% keyword in the plot command (added in the example above (Jan)) - or by setting it globally with:
<code idl>
  !X.MARGIN = [10,3]
</code>
X margin gives the width of the margin in characters on the left (XMARGIN[0]) and right (XMARGIN[1]) side of the plot axis.
Similarily the bottom and top margins can be set with
<code idl>
  !Y.MARGIN = [4,2]
</code>
or with the YMARGIN keyword in the plot command.


==== multiple plots ====

To get multiple sub plots arranged in a matrix use !P.multi :

  !P.multi = [ $
    N_remaining , $ ; number of plots remaining on page , a zero starts a new page
    N_column , $ ; number of columns
    N_row , $ ; number of rows
    N_z_stack , $ ; number of plots stacked in z dimension
    top_2_bottom $ ; Flag: first fill columns before proceeding to next, default is to first fill rows ('column major')
     ] 
 
Usage:
  !P.multi = [ 0 , 2, 3 ]
plots come in a new (erased) 2x3 array, next plot will go into top left position, following one in top row second column, etc.
  !P.multi = [ 0 , 5, 5, 1, 1 ]
plots come in a 5x5 array, next plot will go into top left position, next one in second row, first column, etc.
  !P.multi = [ 0, 3, 3 ]
  ...
  !P.multi = [ 5 ]
Next plot goes into row 2 colomn 2




==== Textausgabe in Grafik ====

  XYOUTS, X, Y, String, alignment=alignment, charsize=charsize, width=width, clip=clip, noclip=noclip, ...

print string at pos x,y \\
  * x,y in data coordinates (axes) - but that can be changed with flags /normalized or /device 
    * omitting x,y position will be last (string?) position. in a new plot this is at 0,0 normalized coordinates (i.e. lower left corner)
  * alignment is a floating point number: 
    * 0.0 : left bounded, i.e. lower left corner of text is at x,y or text starts at x,y
    * 0.5 : centered, i.e. center of text at baseline is at x,y
    * 1.0 : right bounded, i.e. lower right corner  is at x,y
    * yon can also give numbers below 0 or above one - if you know what to do with that
  * charsize is size of character in multiples of character base size (independent of !P.charsize)
  * charthick is linethicknes of drawn fonts (Hershey...)
  * width returns the width of the plotted text in normalized coordinate units, if no other parameter is given the text is not printed.
  * clip=[x0,y0,x1,y1] defines a rectangle where to clip the output text. You need to set noclip=0 to take effect
  * noclip a flag that specifies whether the current clipping region (-> parameter clip or !P.clip) is used, default is noclip=1 i.e. you have to give noclip=0 to take clipping effect.
  * more details [[https://www.harrisgeospatial.com/docs/XYOUTS_Procedure.html|xyouts@harris]]

  * it might be useful to give positions in units of character height or width they can be determined from the system in device coordinates:

  ; char height and (average) width in device coordinates 
  ch = float(!D.y_ch_size) 
  cw = float(!D.x_ch_size) 
  ; ... in normalized coordinates  
  ch = float(!D.y_ch_size) / !D.y_size
  cw = float(!D.x_ch_size) / !D.x_size
  ; ... in data coordinates (only in linear coordiante system !)
  ch = float(!D.y_ch_size) / (!D.y_size * !Y.S[1])
  cw = float(!D.y_ch_size) / (!D.y_size * !X.S[1])

  * it seems as if for pixel based devices it is %%!D.y_ch_size = 12%% pixel. To achieve a certain charsize q_ch relative to bitmap size do:

  ; charsize in normalized units
  ; ch = !P.charsize * float(!D.y_ch_size) / !D.y_size
  ; I want ch = q_ch% <=> q_ch = charsize * float(!D.y_ch_size) / !D.y_size * 100
  ; => charsize = q/100. * !D.y_size / float(!D.y_ch_size)
  q_ch = 3. ; percent of bitmap height
  !P.charsize = q_ch/100. * !D.y_size / float(!D.y_ch_size)

==== Symbol in Grafik ==== 

Um ein einzelnes Symbol an bestimmter Stelle in einem Plot zu positionieren:

  PLOTS,  xpos, ypos , psym=<symbol id>, ... 

==== Fonts ====

there are three font systems in IDL:
  * !P.Font = -1 [[https://www.nv5geospatialsoftware.com/docs/Using_Hershey_Vector_Fon.html|Hershey-Vector Fonts]] (scalable and rotatable Vectorfonts, default of idl)
  * !P.Font = 0 [[https://www.nv5geospatialsoftware.com/docs/using_device_fonts.html|'device Fonts']] - in the Postscript device these are Postscript Fonts - they are better readable than true-type fonts
  * !P.Font = 1 [[https://www.nv5geospatialsoftware.com/docs/using_truetype_fonts.html|True Type Fonts]] (TTF if installed). With TTF some of the below described formatting stuff is not possible - I (jan) find them rather 'pixelig'


=== Select font ===

IDL provids the possibility to change the font within the %%xyouts%% routine or the %%[x,y,z]title%% of plots.

  * put an escape sequence of  the form %%!<command>%% in the string to be printed
  * different fonts can be selected with  %%!<digit>%% here comes a selection
    * !3 is the Hershey Helvetica Font which is the default. IDL calls it "simplex Roman" - obviously the people at IDL mixed up the names
    * !6 is the Hershey Font "complex Roman" where every lines is made 3times to make the letters more bolder.
    * !4 is the Hershey Font "simplex Greek" greek letters single line: greek letters in principle in standard encoding, i.e. a ist alpha etc. but q instead of r is rho ... there is somewhere a shift in the encoding table! (see [[https://www.nv5geospatialsoftware.com/docs/Using_Hershey_Vector_Fon.html|"Using vector fonts"]] ).
    * !7 is Hershey Font "complex Greek" greek letters double lines  
    * !9 oder !M is "special Math" - Integral and other special Math-signs this font is used for only one letter
    * !X returns to the 'Entry font'
    * for examples of the available fonts see IDL help under [[https://www.nv5geospatialsoftware.com/docs/Using_Hershey_Vector_Fon.html|Idl help "Using vector fonts"]]. There you will also see the encoding. You can use the indices from this tables directly. To e.g. print the degree sign which you find in the default font 3 under octal '260' use:\\ %%  xyouts, x,y, '15'+string("260b)+'C'%%\\ (the notation is correct: quotation marks do not close and b stands for ocal !)((IDL has some odities he !))

Examples with simplex greek
  * a small greek phi\\ %%!4f!X%% (take font greek simplex (!4) set an f return to default font (!X))
  * a large greek phi\\ %%!4F!X%% 
You may also use the complex greek (2lines per stroke to make letter more bold)
  * a greek psi\\ %%!7w!X%% 

One may define constants for the different letters:
  c_greek_font = '!4' ; simplex greek
  ;c_greek_font = '!7' ; complex greek
  c_alpha   = c_greek_font+'a!X'
  c_beta    = c_greek_font+'b!X'
  c_gamma   = c_greek_font+'c!X'
  c_delta   = c_greek_font+'d!X'
  
  c_epsilon = c_greek_font+'e!X'
  c_zeta    = c_greek_font+'f!X'
  c_eta     = c_greek_font+'g!X'
  c_theta   = c_greek_font+'h!X'
  
  c_iota    = c_greek_font+'i!X'
  c_kappa   = c_greek_font+'j!X'
  c_lambda  = c_greek_font+'k!X'
  c_mu      = c_greek_font+'l!X'
  
  c_nu      = c_greek_font+'m!X'
  c_xi      = c_greek_font+'n!X'
  c_omikron = c_greek_font+'o!X'
  c_pi      = c_greek_font+'p!X'
  
  c_rho     = c_greek_font+'q!X'
  c_sigma   = c_greek_font+'r!X'
  c_tau     = c_greek_font+'s!X'
  c_upsilon = c_greek_font+'t!X'
  
  c_phi     = c_greek_font+'u!X'
  c_xi      = c_greek_font+'v!X'
  c_psi     = c_greek_font+'w!X'
  c_omega   = c_greek_font+'x!X'
  
  c_infinity= c_greek_font+'y!X'
  
  ; usage
  xyouts, 0,0, c_phi+'('+ c_theta+') = '+ c_delta


=== Hershey-Vector Fonts ===

Hersehy vector fonts are one of the oldest computer fonts. They were developed in the late sixties by A.V. Hershey, see [[http://en.wikipedia.org/wiki/Hershey_font|wikipedia]], [[http://paulbourke.net/dataformats/hershey/| Paul Bourke]] or [[https://www.nv5geospatialsoftware.com/docs/Using_Hershey_Vector_Fon.html|Idl help "Using vector fonts"]]. They were designed for a plotter i.e. a device which draw lines with a pen. The letters consist solely of polygons or straight lines. To achieve bold(er) letters two or three lines are drawn close to each other. The respective fonts are called then duplex or triplex respectively. For examples see [[https://www.nv5geospatialsoftware.com/docs/Using_Hershey_Vector_Fon.html]]. To find out the byte code of the symbols shown in the tables add the (decimal) values at the left and top of the respcetive rows and columns. To get e.g. an $\alpha$ , i.e. the symbol in the row labelled with 96 and the column labelled with 1 you can use %% string(byte(96+1))%%.

==== postscript fonts ====

setting the post script device up with

  SET_PLOT,'ps' ; use postscript device
  !P.font = 0 ; use postscript device fonts 
  DEVICE, /ENCAPSULATED, /portrait, /COLOR, /HELVETICA , filename='test_font.eps'
  
  ; set some font indices (overwrite the original ones )
  ; index must be in [3..19]
  ; odd indices are normal , even are italic/oblique
  ; we need one without serifs (helvetica) one with (times) a typerwriter (courier) and symbol 
  DEVICE, /HELVETICA ,           FONT_INDEX=3
  DEVICE, /HELVETICA , /OBLIQUE, FONT_INDEX=4
  DEVICE, /TIMES ,               FONT_INDEX=5
  DEVICE, /TIMES ,     /ITALIC,  FONT_INDEX=6
  DEVICE, /COURIER ,             FONT_INDEX=7
  DEVICE, /COURIER ,   /OBLIQUE, FONT_INDEX=8
  DEVICE, /SYMBOL ,              FONT_INDEX=9 
  
  ; char height in normalized coordinates
  ch = float(!D.y_ch_size) / !D.y_size
  
  ; some letters...  
  abc_rst = 'abcdefghijklmnopqrst...'
  ; character size factor
  cs = 2.5
  ; sample output: switch to font #4 with !4 and back to 'standard' with !X
  ; gor through fonts 3-9 defined above
  for i = 3, 9 do begin
    si = string(i,f='(I0)')
    s = '!'+si+'font #'+si+':'+abc_rst+'!X'
    print, s
    xyouts, 0.01, 1-(i-3+1)*cs*ch, s, charsize=cs, /normal
  endfor
  
  device, /close


see the end of the [[https://www.nv5geospatialsoftware.com/docs/Using_Device_Fonts.html |  Device fonts ]] section.
           

==== type setting ====

You may generate super- and sub-scripts and further primitve type setting command based on 
[[http://chu.berkeley.edu/dokuwiki/chu:soft:fontart|Grandle and Nystrom (1980)]]((Grandle and Nystrom, 1980, A method for usage of a large calligraphy set under RSX-11M. Grandle, R.E. and Nystrom, P.A. Proceedings of the Digital Equipment Computer Users Society, November 1980, 391-395.)). This is also used in other software packages see e.g. the documentation of  <del>[[http://www.nag.co.uk/visual/IE/iecbb/DOC/html/unix-ref/stan/nagtext.htm|NAG]] e.g.  [[http://www.ualberta.ca/AICT/RESEARCH/NAG/Explorer5doc/html/unix-ref/stan/nagtext.htm|here]] or [[http://www.nag.co.uk/visual/IE/iecbb/DOC/html/unix-ref/stan/nagtext.htm|here]]</del> --- sorry all the links before are dead. 
Here are some still alive NAG docs pages:
[[http://bh0.phas.ubc.ca/Doc/explorer/doc/html/doc/ref/stan/nagtext.htm|NagText]], 
[[https://sites.cc.gatech.edu/scivis/iris/doc/ref/stan/nagtext.htm|NAGText@gatech.edu]], 
[[https://sites.cc.gatech.edu/scivis/iris/doc/ref/stan/naggraph.htm|NAGGraph@gatech.edu]].

The following commands are available (see also IDL help [[https://www.nv5geospatialsoftware.com/docs/Embedded_Formatting_Comm.html|Embedded Formatting Commands]])((find it under 'fonts' in the help)):

  * supersuperscript  !U ('Up' faktor 0.62)
  * subscript    !D ('Down' faktor 0.62)
  * subsubscript !L ('Low' faktor 0.62)
  * index        !I (faktor 0.44)
  * exponent     !E (faktor 0.44)
  * back to Normal !N
  * numerator position !A ('above' division line)
  * denominator position !B ('below' division line)
  * new line !C ('carriage return')
  * special characters !Z(b1[,b2[,...]]) where b1,b2,... are the 16 bit hexadecimal unicode codes of the desired character[s]
  * Save position !S
  * Restore Position !R 
  * exclamation mark !!

Most of these commands only affect the vertical position except for %%!R%% which restrores also the horizontal position.

Examples:

  * **super and subscripts**
    * $A^b$ (A to the power of b) \\ %%A!Eb!N%% \\(A !Exponent b !Normal ) 
    * To get both super and sub script directly after your letter you need to save the horizontqal position: \\ $A^b_i$   (A superscript b index i)\\ %%A!S!Eb!N!R!Ii!N%% \\ (A !Save !Exponent b !Normal !Restore !Index !Normal)\\ If super and subscript have different lengths you first need to set the shorter one.

  * **X-axis Annotation** in two lines: use code !C for a newline
<code>
  plot, time, y, $
    xtitle = '!CDate (!S!Add.mm!R!B YYYY !N)', $
    xtickformat='(C(CDI02,".",CMoI02,"!C",CYI04))', $
    ...
</code>
    * the xtitle starts with a newline (!C) because IDL does not recognize that the annotations have two line now.\\ !S!A saves the position after the "(" and writes "dd.mm" as nominator,\\ !R!A returns to the position after "(" and writes the denoniator \\ !N ends the nominator and returns to normal before the closing ")" is set. \\ This generates something like:
<code>
  ------+-----------+----------> 
      24.12       24.01
      2019        2020
</code>
    * something similar can be achieved by using  %%  xtickformat='(C("!S!A",CDI02,".",CMoI02,"!R!B",CYI04,"!N"))' ...%% with no centering of the YYYY content

  * **special characters**: 
    * the **degree sign**: 
      * a small letter 'o' in the exponent: $15^oC$ : %%15!Eo!NC%%
      * degrees as the special character with octal code 260: %%string("260b)%%
      * via its decimal code 176:  %%string(byte(176))%%
      * or via its hexcode 0xB0:  %%!Z(B0)%%

    * special characters: greek letters in simplex greek and complex greek
      * $\beta$: %%!4b!X%% (switch to simplex greek (!4), draw b, go back to previous font (!X))
      * $\beta$: %%!7b!X%% (switch to complex greek (!7), draw b, go back to previous font (!X)) \\ note that encoding of the greek letters follows only for the four first letters the latin alfabet i.e. a hives $\alpha$ b give $\beta$, but r gives not $\rho$ you need q instead, and s not $\sigma$, you need t instead.

    * cdot (center dot for multiplication)
      * $\cdot$: %%'!9'+string(byte('056'O))+'!X'%% \\ select special and math font (!9), select the cdot symbol (string(byte(...))), return to previous font (!X)

    * integral sign with integration boundaries
      * $\int_a^b dx$ : %%'!9'+string(byte(105))+'!X!S!Da!N!R!Ub!N dx'%% \\ select special math font (!9), integral sign (string(byte(105))), restore previous font (!X), save position (!S), lower boundary (!Da!N), restore position, upper boundary (!Ub!N), dx. \\ the inteegral sign used here is rahter small alternatively one can use string(byte(73)) also in the MAth and special font (!9) but it is rather large ...

    * a vector and its absolutum
      * $\vec{v}$: %%  '!Sv!R!9!U'+string(byte(54))+'!N!X'%% \\ (save position(!S), draw letter v, go back to start position (!R), select special font (!9), upper position (!U), arrow (string(byte(54))), normal position (!N), previos font (!X) )
      * $||\vec{v}||$:  %%  '!9'+string(byte(35))+'!X'!Sv!R!9!U'+string(byte(54))+'!N'+string(byte(35))+'!X'%% \\ where the aboslutum symbol is generated with %%'!9'+string(byte(35))%% 

  * more complex constructions show the limits: 
    * **fraction**  $\frac{a}{1-x}$  \\  %%!S!A a !R!B1-x!N%% \\ i.e.: save position (!S), draw nominator (!A), go back to saved position and start denominator (!R!B), return to normal when done (!N at end)  \\ as there is no automatic centering the %%a%% in the nominator is surrounded here with two spaces which is not exact. And there is no fraction bar: you have to do it on your own which means you need to know the length of the object. This could look like this (with no guarantee):
<code>
  xyouts, x,y , '!S!A a !R!B1-x!N' , width=wdt 
  plots , [x, x+wdt ] , [ y+dy , y+dy ]%%
</code>
    * The position of the fraction line must be corrected a bit upward therefore the %%dy%%((y is the position of the base line on which printing occurs. The fraction line is a bit less then half a character height above)). You may try to calculate %%dy%% but at a certain point it becomes incredible complex ...

For more examples have a look at [[https://www.nv5geospatialsoftware.com/docs/formatting_command_examp.html#fonts_2203613354_1050540|Formatting Command Examples]].



==== graphics -> screen ====

With a call to plot (and similar functions) a grafic window is opened (except you have changed the output device, see below).
This x-terminal (linux/linux) does not perform a refresh when covered by another window and it uses only 256 colors.

For an automatic refresh set the retain flag to two:

  SET_PLOT,'X'
  DEVICE, RETAIN=2

This works only if the x-window is opened for the first time. If a x-terminal is already open close it before you call the device.

If you want "true colors" (i.e. 24bit colors) use:

  SET_PLOT,'X'
  DEVICE, /DECOMPOSED

Colors are then defined by hexadecimal numbers. Red, green, blue and combinations can be defined by:

  red   = '0000ff'x
  green = '00ff00'x
  blue  = 'ff0000'x
  yellow= '00ffff'x
  magenta='ffff00'x
  cyan  = 'ff00ff'x
  black = '000000'x
  white = 'ffffff'x
  grey  = '888888'x



==== Farbe, colortables ====

Farbe wird in idl üblicherweise als Index in eine Farbpallete (=colortable) interpretiert. Diese Tabellen können maximal 256 Einträge enthalten. Einige plotdevices können aber auch true color (also 3*8Byte=24bit farbwerte). 

Alle Ausgaben in das Plotdevice können farbig erfolgen - dies geschieht mit der %%color = <wert>%% option, wobei %%<wert>%% entweder ein index in eine Farbpalette ist oder ein true colorwert.


=== colortables ===

From version 8 on IDL takes truecolor rgb values which you give as 3byte integer values. 
You can write them as Hex-numbers like color='0000ff'x, color='00ff00'x and color='ff0000'x for red, green and blue, respectively. 

Before version 8 IDL used colortables, the color appearing on the screen was defined in a table or palette, i.e. a list with 256 entires for the r,g,b truecolor values. There are several predefined colortables (currently 73). You can load them with
<code idl>
LOADCT, ct_index
</code>
where ct_index is the index of the color table. 
A command %%PLOT,x,y,color=4%% plots data with the color defined in the palette for index 4. 
The default palette ist a grayscale grayscale with 256 entries.

To print a list of the color table names use
<code>
  loadct, get_names=s
  for i = 0, n_elements(s)-1 do print, i, ':', s[i]
</code>
Current list is 
       0:B-W LINEAR
       1:BLUE/WHITE
       2:GRN-RED-BLU-WHT
       3:RED TEMPERATURE
       4:BLUE/GREEN/RED/YELLOW
       5:STD GAMMA-II
       6:PRISM
       7:RED-PURPLE
       8:GREEN/WHITE LINEAR
       9:GRN/WHT EXPONENTIAL
      10:GREEN-PINK
      11:BLUE-RED
      12:16 LEVEL
      13:RAINBOW
      14:STEPS
      15:STERN SPECIAL
      16:Haze
      17:Blue - Pastel - Red
      18:Pastels
      19:Hue Sat Lightness 1
      20:Hue Sat Lightness 2
      21:Hue Sat Value 1
      22:Hue Sat Value 2
      23:Purple-Red + Stripes
      24:Beach
      25:Mac Style
      26:Eos A
      27:Eos B
      28:Hardcandy
      29:Nature
      30:Ocean
      31:Peppermint
      32:Plasma
      33:Blue-Red
      34:Rainbow
      35:Blue Waves
      36:Volcano
      37:Waves
      38:Rainbow18
      39:Rainbow + white
      40:Rainbow + black
      41:CB-Accent
      42:CB-Dark2
      43:CB-Paired
      44:CB-Pastel1
      45:CB-Pastel2
      46:CB-Set1
      47:CB-Set2
      48:CB-Set3
      49:CB-Blues
      50:CB-BuGn
      51:CB-BuPu
      52:CB-GnBu
      53:CB-Greens
      54:CB-Greys
      55:CB-Oranges
      56:CB-OrRd
      57:CB-PuBu
      58:CB-PuBuGn
      59:CB-PuRd
      60:CB-Purples
      61:CB-RdPu
      62:CB-Reds
      63:CB-YlGn
      64:CB-YlGnBu
      65:CB-YlOrBr
      66:CB-BrBG
      67:CB-PiYG
      68:CB-PRGn
      69:CB-PuOr
      70:CB-RdBu
      71:CB-RdGy
      72:CB-RdYlBu
      73:CB-RdYlGn
      74:CB-Spectral


Die aktuell geladene Farbpalette kann mit 
  TVLCT, rgb , /Get
in das array %%rgb%% geladen werden, oder mit 
  TVLCT, red, grn, blu , /Get
in die arrays red, grn, und blu. 

Der folgende Beispielcode liest die aktuelle Farbpalette und definiert ein paar Farben um und macht diese **neue Palette** zur aktuellen:

<code idl>
; get actual colortable
TVLCT, rgb , /Get
; size of this colortable
n_col = n_elements( rgb[*,0] )
; define some indices into the table
col_black = 0
col_red = 1
col_green = 2
col_blue = 3
col_yellow = 4
col_gray = 5
col_white = n_col-1
; define the rgb components of the colors related to this indices
rgb[ col_red   , * ] = [   0 ,   0 , 255 ]
rgb[ col_green , * ] = [   0 , 136 ,   0 ] 
rgb[ col_blue  , * ] = [ 255 ,   0 ,   0 ]
rgb[ col_yellow, * ] = [   0 , 255 , 255 ]
rgb[ col_gray  , * ] = [ 136 , 136 , 136 ]
; IDL uses the first and the last entry of the colortable 
; as background and primary drawing color i.e. as black an white 
; redefine it here ...
rgb[ col_black , * ] = [   0 ,   0 ,   0 ] 
rgb[ col_white , * ] = [ 255 , 255 , 255 ]
; make this new colortable active
TVLCT, rgb
...
; plot red axis with red line
plot, x,y , color=col_red
...
</code>

Ähnlich wirkt der IDL befehl [[http://www.physics.nyu.edu/grierlab/idl_html_help/T7.html#wp137842|TEK_COLOR]] der allerdings keine Namen definiert. First nine colors are: 0=black, 1=white, 2=red, 3=green, 4=blue, 5=cyan, 6=magenta, 8=orange.

Klassiche Farbpalleten die man auf diese Weise definieren könnte sind z.B. die [[http://en.wikipedia.org/wiki/Color_Graphics_Adapter#Color_palette|CGA-Palette]], die [[http://en.wikipedia.org/wiki/Enhanced_Graphics_Adapter#The_EGA_color_palette|EGA-palette]]  oder die [[http://en.wikipedia.org/wiki/Video_Graphics_Array#The_VGA_color_palette|VGA-pallette]].

Einen Farbverlauf kann man ähnlich definieren. Der folgende code definiert ein **Spektrum** in dem, orientiert an vereinfachten Empfindlichkeitskurven der Farbrezeptoren des Auges die Kurven für rot, grün und blau gesetzt werden:

<code idl>
; calc. a spectrum starting with red going via yellow and green to blue and violet
; get current color table
TVLCT, RED, GRN, BLU , /Get
; size of this colortable
n_colors = n_elements( RED )
; maximum value of one channel
max_rgb = 255
; the very first and last points (col=0 and col=N-1) are excluded because idl uses them as background and foreground colors
for col = 1 , N_colors-2 do begin
  RED(col) = fix((cos(2.*!pi*0.8*(float(col-1)/(N_colors-3)    ))+1.0)/2.0*max_rgb)
  GRN(col) = fix((cos(2.*!pi*    (float(col-1)/(N_colors-3)-0.3))+1.0)/2.0*max_rgb)
  BLU(col) = fix((cos(2.*!pi*0.5*(float(col-1)/(N_colors-3)-1.0))+1.0)/2.0*max_rgb)
  endfor
; load this color table
TVLCT, RED, GRN, BLU

</code>


=== true color, direct color ===

In the **'true color'** or **'direct color'** system are the 3 lowest bytes of an integer number interpreted as the intensity values for red, green and blue. If you write this as a hex number you may define rather easily arbitrary colors. Here are some basic colors as true colors (compare also with 'html colors' eg in [[http://www.w3schools.com/html/html_colors.asp]]):
<code idl>
col_black  = '000000'x
col_red    = '0000ff'x
col_green  = '008800'x
col_blue   = 'ff0000'x
col_yellow = '00ffff'x
col_gray   = '888888'x
col_white  = 'ffffff'x
</code>

To use true colors the plotdevice must be configured accordingly. This is done with the /decomposed keyword, e.g. for UNIX X-windows device:
<code idl>
SET_PLOT, 'X'
DEVICE, /decomposed
</code>

One may also define colors directly as 24bit numbers from red, green and blue values. 
Assume that red, grn and blu are byte values (8-bit numbers) in the range 0..255, the respective '24bit true color' values is then:
<code idl>
true_color_rgb = long(red) + ishft(long(grn),8) + ishft(long(blu),16)
</code>

or if you have your red, grn blu values as floats in the range 0..1
<code idl>
true_color_rgb = long(red*255) + ishft(long(grn*255),8) + ishft(long(blu*255),16)
</code>




==== Plot system variables ====

!P is a structure containing system parameters describing the **plot**

!P.MULTI -> **multiplot** 
  * !P.MULTI[0]  number of plots remaining on the page
  * !P.MULTI[1]  plot columns per page.
  * !P.MULTI[2]  rows of plots per page. 
  * !P.MULTI[3]  plots stacked in the Z dimension.
  * !P.MULTI[4]  is 0 to make plots from left to right (column major), and top to bottom,\\ 1 to make plots from top to bottom, left to right (row major).


!P.charsize -> relative size of all characters, default: 1=normal size
  * 2 would double the size of all text in the plot
  * 0.75 would reduce the size of all characters by about 25%

!D is a structure containing parameters about the current plot-**Device**\\
!D.x_size and !D.y_size give the width of the device in device coordinates\\
!D.x_ch_size and !D.y_ch_size give the base size of characters in device coordinates\\
=> !D.y_ch_size / !D.y_size gives the base size of characters in normalized coordinates\\



!X, !Y and !Z are structures containing parameters determining the appearance of the X,Y, and Z **axes**\\
!X.Margin = [left, right] give the width of the left and right margin between the axes and the plot region in characterwidths\\
!Y.Margin = [bottom, top] give the width of the left and right margin between the axes and the plot region in characterheights\\

Some variables contain information about the **last plot**

!X.S and !Y.S give the scaling factor for conversion from data to normalized coordinates. It is
  x_norm = !X.S[0] + !X.S[1] * x_data
  x_norm = !X.S[0] + !X.S[1] * log_10(x_data)
for linear and logarithmic scaling respectivley

!X.crange = fltarr(2) contains the range used by the x-axis.

!X.window = fltarr(2) contains the normalized x-coordinates of start and endpoint of the axis.

!X.region = fltarr(2) contains the normalized x-coordinates of the plotting region

Every single plot goes into a **region** which should contain the plot itself, the axis labels, the axis titles and the plot title. !X.Margin determines the distance of the x-axis to the bottom margin of this plot region.

If you want to know where ticks are made then use keywords xtick_get and ytick_get when calling plot:
  plot, x, y, xtick_get=xticks, ytick_get=yticks 

See also the [[http://www.harrisgeospatial.com/docs/Graphics_System_Variable.html|Graphics System variables]] in the IDL help.
A (horrible) description of the different coordinate systems of IDL can be found under [[http://www.harrisgeospatial.com/docs/Coordinate_Conversions.html| Direct Graphics Coordinate Systems ]].


==== graphics function ====

For the pixel based Windows, X and Z Device it is possible to define a bit-operation for every pixel drawn. I.e the color bits for every pixel of a line are first connected with the background of its position. This can be used to make lines contrasting with the underlying background.
The operation can be defined with the SET_GRAPHICS_FUNCTION parameter of the DEVICE function :

  DEVICE, GET_GRAPHICS_FUNCTION = old_gfct, SET_GRAPHICS_FUNCTION = 6

with graphics_function = 6 = GXxor an XOR operation between pixel color and the content of the pixel. The GET_GRAPHICS_FUNCTION variable stores the current state.
The default is GXcopy = 3 - i.e 
A list of the possible grafics functions can be found under idl help -> [[https://www.nv5geospatialsoftware.com/docs/DEVICE_Procedure.html#devices_517620971_499029| Keywords Accepted by the IDL Devices -> SET_GRAPHICS_FUNCTION ]]

The result of this pixel operations can give weird results ...

==== transparency ====

IDL does not know transparent overlays but on pixel based devices you can implement it with TVRD and TV.
Lets assume you want to put a bitmap BMP with size Nx, Ny at a position x0 y0 (pixels) with a tranparency factor of q_tra:

  ; get current content of canvas as a bytarr( N_channel, Nx, Nx )
  bgnd = TVRD( x0,y0, Nx, Ny, true=1 )
  ; blurr it a bit
  N_blurr = 8
  bgnd = SMOOTH( bgnd, [ 1, N_blurr, N_blurr ] , /edge_truncate )
  ; create a new bitmap by combining bmp and bgnd 
  bmp = byte( q_tra * float(bmp)  +  (1-q_tra) * float(bgnd) )
  ; put the new bmp into the canvas
  TV, bmp, x0, y0, true=1, /device

This is implemented in mk_legend.pro 
==== coordinate transform ====

in classical 2D plots IDL knows three coordinate systems. In plot commands as plot, oplot, xyouts, plots you can refer to them via key words
  * **/data**    the coordinate system defined by the current axes
  * **/normal**  normalized coordinates within the axes running from zero to one
  * **/device**  device coordinates in pixel or postscript points

to convert between them idl uses several parameters (see [[http://www.physics.nyu.edu/grierlab/idl_html_help/plotting5.html#wp952064|Direct Graphics Coordinate Systems ]]):
  * **/data to /normal** is determined via !X.S :  \\ x_norm = !X.S[0] + !X.S[1]*x_data for !X.type=0 (linear) and \\ x_norm = !X.S[0] + !X.S[1]*log(x_data) for !X.type=1 (logarithmic)  \\ and the respective forms for y (-> !Y.S) and z (-> !Z.S)
  * **/normal to device** is determined via !D.X_VSIZE and !D.Y_VSIZE: \\ x_device = x_norm*!D.X_VSIZE




==== Postscript ====

If you want scalable high quality plots with great fonts use encapsulated postscript (eps).
Postscript is a vectorgrafics format and encapsulated means that the grafic can be easily included in other documents (if your software can deal with eps :-)).

But be aware that if your plot is made up by many plot commands the result will be an incredible large file. Nevertheless you can convert eps to pdf (which is still a vector grafics format but has an internal compression, [[internal:administration:idl&#postscript_-_pdf|see below]]) or to png ([[internal:administration:idl&#postscript_-_png|see below]]). If an eps files is very large the pdf will be large (about 1/4) and a png will be 'normal' (what you may used to). //Never use jpeg// for line plots or everyting which is not a photograph. 


to produce an **encapsulated postscript** plot:
  SET_PLOT,'ps' ; select postscript device
  !P.font = 0 ; select postscript fonts
  DEVICE, /ENCAPSULATED, /portrait, /HELVETICA, /COLOR ; basic setup for colored EPS 
  DEVICE, /COLOR, /decomposed ; in idl 8 is true color = /decompose possible 
  ; DEVICE, FILENAME="idl_plot.eps" ; put a filename here if you do not ike the default from idl ('idl.ps') but note that this name is reset with every DEVICE,/CLOSE call
  ; DEVICE, XSIZE=17.8, YSIZE=12.7 ; size in centimeter ! these are the default values corresponding to a 640*480 pixel image printed at 96pixel per inch...
  DEVICE, XSIZE=29.4, YSIZE=21.0 ; size in centimeter ! DIN A4
  DEVICE, /SYMBOL ,  FONT_INDEX=4 ; make greek letters available as font #4
  
  ; device keyword output is intended to add own postscript code,
  ; BUT it generates also ps-code for a new page.
  ; See also comment in http://www.physics.wisc.edu/~craigm/idl/down/pxperfect.pro
  ; DEVICE, OUTPUT = '%%Title: plot by Jan Schween'
  
  ; in post script we need thicker axis lines
  !X.thick = 3
  !Y.thick = 3
  
  ; plot in black on white background
  !P.color = '000000'x
  !P.background = 'ffffff'x
  
  ; for a good readable plot charsize should be larger
  !P.charsize = 2.
  
  ...
  plot commands
  ...
  
  DEVICE,/CLOSE
  
  ; if we did not provide a plotname with DEVICE, plotname+'.eps' we can rename it here...
  SPAWN, 'mv idl.ps '+plotname+'.eps' ; rename outputfile


  * mit /portrait wird das koordinatensystem nicht gedreht und die **bbox stimmt**...\\
  * mit /landscape wird das koordinatensystem gedreht und die **bbox stimmt nicht mehr** ...\\

in idl 6.4 ist unter postscript farbe nur aus einer farpalette mit 256 enträgen möglich. Mit:
  SET_PLOT,'ps'
  DEVICE, ... /COLOR, BITS_PER_PIXEL=8 ...
ist es möglich 2 dimensionale arrays intarr(Nx,Ny) als farbige bitmaps in PS darzustellen.
Die Farbe ergibt sich dabei aus der aktuellen colortable - die hat maximal 256 einträge also einen 8Bit langen indexraum.

Im Postscript device kann man zwar true color jpegs darstellen:
  READ_JPEG, filename, img
  TV, img, TRUE=1
Allerdings wird das jpg bild anscheinend auf eine 8Bit = 256Farb-Palette reduziert...

In idl 8.2 kann man mit true color arbeiten:
  SET_PLOT,'ps'
  DEVICE, /COLOR, /decomposed, ...
ermöglicht die verwendung von true color - z.B.:
  plot, x,y, color='00ffff'x 
benutzt Gelb (=max red + max green) als Farbe.


**Device coordinates** in the postscript device are in 1/1000cm ([[http://www.physics.nyu.edu/grierlab/idl_html_help/devices13.html#wp451367|1000pixel per cm]]). To make a plus sign at 1cm,1.5cm use

  plots, 1*1000,1.5*1000, psym=1, /device


** Postscript - non encapsulated **:
  SET_PLOT,'ps'
  DEVICE [  ,/portrait | ,/landscape ]
  ...
  plot commands
  ...
  DEVICE,/CLOSE
  
(-> .../idl/eps_test/ps_landscape.pro)\\
  * IDL produces per default a 'seascape' plot erzeugt (in gv muss man seascape oder landscape+swap landscape anwählen)
  * /landscape thus has no meaning !
  * /portrait turns the coo-system by 90° to the left ...
  * in contrast to /ENCAPSULATED the bbox is ok then !
  * for details of the postscript device see [[http://www.physics.nyu.edu/grierlab/idl_html_help/devices13.html#wp451367| here ]]





Basic commands, select device driver 
  SET_PLOT, Device

e.q. for Postscript 
  SET_PLOT,'ps'
or for colored EPS 
  SET_PLOT,'ps',/ENCAPSULATED,/COLOR

with only the SET_PLOT command output will be send to the file **'idl.ps'**\\
this can be changed by
  DEVICE, FILENAME='<filename>', ...

the plot-file is then finalized with
  DEVICE, /CLOSE
or
  DEVICE, /CLOSE_FILE

  


=== Difference between postscript (*.ps) and encapuslated postscript (*.eps) ===

  * postscript is a vector grafics format the predecessor of pdf.
That means grafs can be enlarged or shrinked without loosing quality.
  * there is no big difference between PS and EPS:
    * EPS means that the graph can be imported into other applications. 
    * That means only one single page is allowed - no multiple pages
    * the graf //must have// an information about its size - the bounding box (bbox)


==== postscript -> png ====

Convert postscript to png with [[http://www.imagemagick.org/|imagemagick]]:

  convert -density 120 <ps-filename> <filename>.png

density 120 means 120 dpi, increasing this values increases quality and size of the png.
You must have [[http://www.ghostscript.com/|ghostscript]] installed. 
For more control you may want to directly use the ghostscript call for this:

  gs -dSAFER -dBATCH -dNOPAUSE -sDEVICE=png16m -r120 -dTextAlphaBits=4 -dGraphicsAlphaBits=1 -dEPSCrop -sOutputFile=<filename>.png <ps-filename>

Parameters %%-dSAFER -dBATCH -dNOPAUSE%% makes gs to run this as batch job with no questions to the user,\\
%%-sDEVICE=png16m%% defines true color rbg png as output format, \\
you may use %%png256%% or %%png16%% for 256 or 16 element palette images which give substantilly smaller images but you may experience ditherung artefacts 
(for details on the gs devices see [[https://www.ghostscript.com/doc/current/Devices#PNG|gostscrip webpages]]).
\\
%%-r120%% means 120dpi and gives a 1388x992 pixel images for a DINA4 landscape postscript,  \\
%%-dTextAlphaBits=4 %% ensures smooth characters after rastering, \\
%%-dGraphicsAlphaBits=1%% avoids gridlines after rastering (http://www.idlcoyote.com/ps_tips/psstripes.html), \\
%%-dEPSCrop%% crops the image at the bounding box.

You can execute gs from within IDL with the SPAWN command. As of IDL 9.x you need to prepend a ''unset LD_LIBRARY_PATH'' because IDL try to ensure its own (outdated) libraries for any spawn call:
  SPAWN, 'unset LD_LIBRARY_PATH; gs -dSAFER -dBATCH -dNOPAUSE -sDEVICE=png16m -r120 -dTextAlphaBits=4 -dGraphicsAlphaBits=1 -dEPSCrop -sOutputFile='+plt_name+'.png '+pl_name+'.eps', result, err_result

==== postscript -> PDF ====

One can convert eps files into pdf using the program ps2pdf of the ghostscript package.
The following command creates a pdf file quelle.pdf from the file quelle.eps where the 'papersize'
is adapted to the BoundingBox of the eps file.

  ps2pdf -dEPSCrop quelle.eps

Since ps2pdf is based on ghostscript all [[http://pages.cs.wisc.edu/~ghost/doc/cvs/Use.htm|gs options]] can be used. 

To fit the plot to a standard paper size (i.e. Letter or DINA4 depending on your system settings) use

  ps2pdf -dEPSFitPage quelle.eps


One can also invoke ghost script (gs) to control more precisely what you want.

To create a PDF with a rgb colors use :

  gs -dBATCH -dNOPAUSE -sColorConversionStrategy=RGB -sDEVICE=pdfwrite -dEPSCrop -sOutputFile=quelle.eps.pdf  quelle.eps

To create a PDF based on PDF-A standard (the definitions file PDFA_def.ps should be in the gs/lib dirictory ... no guarantee that this works ...)

  gs -dPDFA=1 -dBATCH -dNOPAUSE -sColorConversionStrategy=RGB -sDEVICE=pdfwrite -sOutputFile=quelle_a.pdf  PDFA_def.ps  quelle.eps


==== SVG ====

Scalable Vector graphics can be displayed by web browsers and might get the standard vector graphic format for web applications

There are converters from eps to SVG:

=== pdf2svg ===

is a simple converter from pdf to svg. Simple because it converts fonts to vector grafics...
The conversion is done in two steps:\\
1: convert eps to pdf
  ps2pdf -dEPSCrop wind_rose_all.eps

2: convert pdf to svg
  pdf2svg wind_rose_all.pdf wind_rose_all.svg

but fonts look really ugly !


=== pstoedit ===

is a package to convert postscript to something. To convert to svg use:

  pstoedit -f plot-svg -psarg "-dEPSCrop" <file.eps>  <file.svg>
  
unfortunately the -dEPScrop option does not work properly. The graph sits in the lower corner and is partially cropped.


==== JPG, PNG and GIF ====

Raster-Oriented (or Bitmap) graphics can be created with:

<code idl>

  ; set the z-device => plot commands go into the Z-Buffer - this is just a bitmap
  SET_PLOT,'Z'
  ; configure this bitmap: true color (24bit=3byte) and size 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=24, /DECOMPOSED, set_resolution=[ 640 , 480 ]
  ...
  plot commands 
  ...
  ; read actual graf into array img as BYTARR( N_color , width , height )
  img = TVRD( true=1 )
  ; write into png image
  WRITE_PNG, 'idl.png', img
  ; or write into jpg image
  WRITE_JPEG, 'idl.jpg', img
</code>

This creates true color graphic files. 

If a 256 color palette shall be used - which reduces the size of the grafics file and which is the default for idl and for gif but also possible with png the code must be changed:

<code idl>
  ; configure the bitmap: 256 color palette (8bit=1byte) and size 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=8, set_resolution=[ 640 , 480 ]
  ...
  ; read actual graf into array img as BYTARR( N_color , width , height )
  img = TVRD( )
  ; read actual palette
  TVLCT, red, grn, blu , /Get
  ; write into png image
  WRITE_PNG, 'idl.png', img, red, grn, blue
</code>

It is also possible to copy an image from the X-Device (i.e. the screen) but this might be disturbed 
by redraw commands or overlap from other windows (see /home/jschween/exp/ideen/idl/mk_png/...) 
=> Use the Z-Buffer instead.


=== animated GIF ===

It is possible to create animated gifs by repeatedly writing to the same gif file with the /MULTIPLE flag set.
GIF knows only palette based color with in maximum 256 colors. We thus either need to work in a 8bit color device, or we quantizie the image colors before providng it to the %%WRITE_GIF%% :

From a 8 bit color device :

<code idl>
  ; plot into virtual device = bitmap
  SET_PLOT, 'Z'
  ; configure this bitmap: 8 bit colors and size is 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=8, set_resolution=[ 640 , 480 ]

  ; get current color palette
  TVLCT, rgb , /get

  ; go through images
  FOR i_img = 0, N_img-1 DO BEGIN

    ; create plot 
    ...

    ; read actual graf into array img
    bmp = TVRD( )
    ; append img to gif-file
    WRITE_GIF, 'idl.gif', bmp,  rgb[*,0],rgb[*,1],rgb[*,2],  /MULTIPLE, delay_time=6, repeat_count=0

  ENDFOR ; i_img

  ; close gif
  WRITE_GIF, 'idl.gif', img, /CLOSE

</code>


From a true color device

<code idl>
  ; plot into virtual device
  SET_PLOT, 'Z'
  ; configure this bitmap: 24 bit colors and size is 640*480 (default)
  DEVICE, SET_PIXEL_DEPTH=24, set_resolution=[ 640 , 480 ]

  ; go through images
  FOR i_img = 0, N_img-1 DO BEGIN

    ; create plot 
    ...

    ; get current bitmap
    bmp = TVRD( true=1 )
    ; quantize colors - derive gif color translation table from first image
    if i_img eq 0 then gif = COLOR_QUAN( bmp, 1, red, grn, blu, get_translation = gif_col_trans ) $
                  else gif = COLOR_QUAN( bmp, 1, red, grn, blu,     translation = gif_col_trans )
    ; append img to gif-file
    WRITE_GIF, 'idl.gif', gif,  red, grn, blu,  /MULTIPLE,  delay_time=6, repeat_count=0

  ENDFOR ; i_img

  ; close gif
  WRITE_GIF, 'idl.gif', /CLOSE

</code>


  * The %%/MULTIPLE%% keyword in the call of %%WRITE_GIF%% creates at the first call a GIF-file but does not close it. During subsequent calls the images are appended. Accordingly the GIF-file must be closed by a final call with the %%/CLOSE%% keyword. 
  * The %%delay_time%% is given in 1/100sec; %%delay_time=6%%/100sec here corresponds to a frame rate of ~17images/sec which is about the frame rate of the old film standard //super8//. Television has usually 25frames/sec corresponding to a delay of 4/100sec.  
  * The parameter %%repeat_count=0%% makes the animation to be repeated forever.
  * Graphics commands must go into a bitmap device in this example the Z-device (it might be possible to do this also in the X-device (UNIX/LINUX), or the win-device (windows). Be aware that the color setup is different than for printer devices (like e.g. postscript): background is black, foreground color (axes of a plot etc.) are white.

=== reading images ===

image files can be read into arrays, displayed, analyzed and manipulated

<code>
  ; read image
  img = read_png( 'img.png' )
  
  ; print out dimensions
  help, img
  
  ; determine dimensions
  img_size = size(img)
  ; number of color channels, columns and rows
  Nc = img_size[1]
  Nx = img_size[2]
  Ny = img_size[3]
  
  ; alternative: query
  qr = query_png( 'img.png' , info )
  Nx = info.dimensions[0]
  Ny = info.dimensions[1]
  
  ; black out lower left corner
  img[ 0 , 0:10, 0:10 ] = 0
  img[ 1 , 0:10, 0:10 ] = 0
  img[ 2 , 0:10, 0:10 ] = 0

  ; display
  x = 10
  y = 10 
  tv, img, x, y , true=1
  
  ; mark center
  plots, Nx/2, Ny/2 , psym=1, /device
</code>



==== back and foreground color ====

The color scheme for bitmap devices as screen (%%SET_PLOT,'X'%% or   %%SET_PLOT,'WIN'%%) or the 'Z'-device (%%SET_PLOT,'Z'%%) is different than for printer devices like e.g. postscript:
  * In bitmaps background is black an drawing is done in white
  * for printers background is white an drawing in black
This can be changed with:
<code idl>
  SET_PLOT, 'Z'
  ; set background color to white (last color in pallete)
  !P.BACKGROUND = 255
  ; set foreground color to black (first color in pallete)
  !P.COLOR = 0
</code>

If you are working with 'true colors' 
<code idl>
  SET_PLOT, 'Z'
  ; set background color to white
  !P.BACKGROUND = 'ffffff'x
  ; set foreground color to black
  !P.COLOR = '000000'x
</code>



==== General remarks about graphic formats ====
  * JPG is good for images with continous color gradients like photographs:
    * (you remember: jp[e]g stands for joint //photographic// experts group).
    * Sharp edges (e.g. lines) in the image create echo like artefacts in jpg.
  * PNG is good for line graphics and images with sharp edges:
    * images with large areas with the same color have good compression rates
    * scaling may lead to moiree efects along edges
    * it is possible to give every pixel in the image a certain transparency
  * EPS / PS ([Encapsulated] Post Script) is good for pure line graphics
    * PS is a vector graphis format. I.E. the file contains plot commands like moveto, lineto, textout, ...
    * if a plot contains a lot of drawing commands the file can become very large
    * bitmap images can be incorporated but are not compressed => files become large...
    * Graphics can be arbitrarily scaled. Scaling is mainly limited by the resolution of the output device.
  * PDF (portable document format) is an extension of PS
    * it has additional compression included
    * EPS and PS images can be converted to PDF ...
    * JPG, PNG etc. can be included in PDF
  * GIF is only good for animations of line graphics - otherwise use PNG or JPG instead !
    * PNG or JPG have more/better colors (GIF has only 256 colors instead of the up to 16million 'true colors' of PNG or JPG)
    * compression is not as good as with PNG or JPG
    * a gif image may have one transparent color i.e. every pixel with this color are totally transparent
    * animation is the only advantage
   

==== polar ====

The plot procedure can work with polar coordinates - you need to set the /polar keyword

  ; define some data points
  N = 100
  theta = findgen(N)/N*360.
  r_max = 7.4 ;-)
  radius = r_max*findgen(N)/N
  
  plot, radius, theta*!DTOR, $
    /polar, $
    xstyle = 1, $
    xrange = [ -r_max, +r_max], $
    ystyle = 1, $
    yrange = [ -r_max, +r_max], $
    /isotropic
  
If you do not fix range and isotropy of the plot it may look weird due to automatic scaling. 

You may want to have axes through the center:

  plot, radius, theta*!DTOR, $
    /polar, $
    xstyle = 1+4, $ ; exact scaling but invisible (box)axis
    xrange = [ -r_max, +r_max], $
    ystyle = 1+4, $ ; exact scaling but invisible (box)axis
    yrange = [ -r_max, +r_max], $
    /isotropic
  
  ; Draw x- and y-axis through 0,0 with ticks to the top or right, respectively.
  AXIS, 0, 0, XAXIS=0
  AXIS, 0, 0, YAXIS=0

You may want to have a radius and angle grid. IDL cant do this for you. You have to do it step by step:

  plot, radius, theta*!DTOR, $
    /polar, $
    xstyle = 1+4, $ ; exact scaling and invisible (box)axis
    xrange = [ -r_max, +r_max], $
    ystyle = 1+4, $
    yrange = [ -r_max, +r_max], $
    /isotropic
  
  ; draw and anotate spokes
  N_spk = 8
  for i = 0, N_spk-1 do begin
    wi = (2*!pi/N_spk)*i
    oplot, /polar, [0,90], replicate( !pi/2-wi, 2 ), linestyle=1
    r = 90.
    xyouts, r*sin(wi), r*cos(wi), string(wi/!dtor,f='(I0)')
  endfor ; i
  
  ; draw and annotate angle circles
  N_w = 181
  w = findgen(N_w)/(N_w-1)*2*!pi
  N_r = 3
  for i = 0, N_r-1 do begin
    ei = (90.*i)/N_r
    oplot, /polar, replicate( 90-ei, N_w), w, linestyle= i gt 0 ? 1 : -1
    xyouts, 0, ei, string(ei,f='(I0)')
  endfor ; i


Some routines like XYOUTS and ARROW do not accept the polar keyword (or ignore it).
You may convert your coordinates to cartesian as above or you can use the CV_COORD function to do this:

  polar_coo = [[angle],[radius]]
  xy_coo = CV_COORD( FROM_POLAR=polar_coo, /DEGREES, /TO_RECT )


==== contour ====

A contour plot with color shading is made with:

  ; prepare number of levels 
  min_data = min( data, /nan )
  max_data = max( data, /nan )
  N_iso = 10
  ; we need one more levels value because level defines interval boundaries, i.e. color[i] goes into interval level[i] ... level[i+1]
  levels = findgen(N_iso+1)/N_iso*(max_data - min_data) + min_data
  ; prepare true colors : e.g. a grayscale from white to black
  gray = long(255*(1-findgen( N_iso )/(N_iso-1)))
  iso_colors = gray or ishft( gray , 8 ) or ishft( gray, 16 ) ; colors as true color 3byte/24bit code
  ; do the contour plot
  contour, data, x_coo, y_coo, $ ; data nad coordinates
    levels = levels, $ ; levels for the colors
    c_color = iso_colors, $ ; 24bit colors
    /cell_fill, $  ; avoid problems with mssing data (nan), etc. (you may use /fill instead but ...)
    xtitle = 'x_coo (units)', $ 
    xstyle = 1, $ ; use exact range of x_coo for axis range
    ytitle = 'y_coo (units)', $
    ystyle = 1 ; use exact range of y_coo for axis range


Routine %%conotur%% takes all the keywords like plot and some more - 
see [[https://www.harrisgeospatial.com/docs/CONTOUR_Procedure.html | help for contour]]
and the descritption of [[https://www.harrisgeospatial.com/docs/graphkeyw.html | graphics keywords]].

Because %%levels%% give the lower border of the interval for the colors we need one more level than colors (%%levels=findgen(N_iso+1)...%%).
For values below level[0] nothing is drawn (no color !), for values above level[N_iso], i.e. the last level the first color is recycled.
To get colors for the values below the lowest level and above the largest level you can add additional colors and levels to what is defined above:

  ; expand coloring levels above and below range(levels)
  iso_colors = [ col_below_min, iso_colors, col_above_max  ]
  levels = [     extrm_lev_min, levels,     extrm_lev_max  ]

you may want to use the most extreme values your system can provide
  extrm_lev_min = -(machar()).xmax
  extrm_lev_max = +(machar()).xmax
  
See function [[https://www.harrisgeospatial.com/docs/machar.html|function machar]].


=== iso lines ===

  contour, data, x_coo, y_coo, levels = levels


If you want iso-lines **with annotations** give the keyword %%/follow%% :

  contour, data, x_coo, y_coo, $
    levels = levels, $
    c_annotation = string( levels , format = '(f5.1)' ), $
    /follow, $
    /overplot

Keyword c_annotations is used here toprovide a format - this is usually not necessary.
You can define also different thicknesses, colors etc. with the other %%c_...%% key words of %%contour%% 
see [[http://www.physics.nyu.edu/grierlab/idl_html_help/C40.html#wp908085| help for contour]].


=== color bar ===

You can use contour also to make a **color bar**:

A **horizontal** color bar at the bottom

  ; charsize in normalized coordinates for position 
  ch = float(!D.y_ch_size) / !D.y_size
  ; the color bar
  contour, [ [levels] , [levels] ]  , $ ; a two row array with just the levels from the contour plot above
    levels, [0,1], $ ; levels on the x-axis and 2 lines in y direction
    levels = levels, $ ; levels for the shading
    c_color = iso_colors, $ ; the colors itself
    /cell_fill, $
    xtitle = 'data', $
    xstyle = 1, $ ; scale x-axis precisely
    xticklen = 5*ch, $ ; ticks must be elongated
    ytitle = '', $ ; no title on y-axis 
    ystyle = 1+4, $ ; do not draw y-axis 
    charsize = 1.0, $ ; smaller charsize ...
    position = [ 0.7 , 2*ch, 0.95, 3*ch ], $ ; postion in lower right corner below axis ...
    /normal, $ ; positon in normalized coordinates
    /noerase ; do not erase beckground - plot over the contour plot

A **vertical** color bar at the right side:

  ; charsize in normalized coordinates for position 
  ch = float(!D.y_ch_size) / !D.y_size
  contour, transpose( [ [levels] , [levels] ] ) , $ ; a two column array with the levels from the contour plot
    [0,1], levels, $ ; two columns in x-direction and levels in y-direction
    levels = levels, $ ; levles for the shading
    c_color = iso_colors, $
    /cell_fill, $
    ytitle = '', $ ; title will be in second y-axis on right side
    ystyle = 1, $
    yticklen = -5*ch, $ ; longer ticks on this axis
    ytickname = replicate(' ', 10), $ ; no annotations on this side 
    xstyle = 1+4, $ ; do not draw x-axis
    charsize = 1.0, $ ; smaller charsize ...
    position = [ 1-4*ch , 7*ch, 1-3*ch, 1-4*ch ], $ ; position at right side ...
    /normal, $ ; position in normalized coordinates
    /noerase ; do not erase background - plot over the contour plot
  
  ; y-axis on the right side 
  axis, yaxis=1, $
    ytitle='data', $
    ystyle = 1, $
    yticklen = -5*ch, $
    yrange=[ min(levels) , max(levels) ], $
    charsize = 1.0, $ ; smaller charsize ...
    /save
    
The position of the vertical color bar in this example is not perfect. You may have to define the position of the contour plot itself to provide space for the color bar.


=== TRIANGULATE ===

In case of data defined on a irregular grid you may need a triangulation, i.e. a kind of ordering of the data in triangles. This triangulation allows you to do CONTOUR plot and further interpolation of this data to a regular grid with TRIGRID.

lets Assume you have data Z at N points defined by coordinates X_coo and Y_coo. To get the triangulation just use:

  TRIANGULATE, X_coo, Y_coo, triangles_idx, boundary_idx
  
triangles_idx is then of type lonarr(3,N_tri) containing the indices to the coordinates of the points defining the N_tri triangles.\\
boundary_idx  is the list of indices of the coordinates of points of the convex boundary.\\
If the boundary of the region covered by your points has concave stretches TRIANGULATE produces triangles which cover them.\\
I found no parameter or function to avoid this or remove them.\\
If your coordinates are given in different, non isotropic units you should **make them isotropic** (e.g. X_coo as geographic coordinates like degrees east, and y_coo in meter above see level -> use %%TRIANGULATE, X_coo*deg_2_meter, Y_coo, triangles_idx%% )

The output can be fed into CONTOUR with

  CONTOUR, Z, X_coo, Y_coo, triangulation=triangles_idx, ...


=== TRIGRID ===
  
You can interpolate triangulated data to a regular grid with TRIGRID:

If you just want a grid with Nx x Ny points

  Z_grid = TRIGRID(  x_coo, y_coo, Z, triangles_idx, Nx=Nx, Ny=Ny, xgrid = x_grid, ygrid = y_grid, missing=!Values.F_NAN )
  
Z_grid is then fltarr(Nx,Ny) and y_grid and y_grid are of type fltarr(Nx) and fltarr(Ny).
  
If you want a section defined by its limits and a step size

  grid_spacing = [ dx, dy ] 
  grid_limits = [ x_min, y_min, x_max, y_max ]
  Z_grid = TRIGRID(  x_coo, y_coo, Z, triangles_idx, grid_spacing, grid_limits, xgrid = x_grid, ygrid = y_grid, missing=!Values.F_NAN )

You can also specify coordinates x_nodes and y_nodes to which to interpolate:

  Z_grid = TRIGRID(  x_coo, y_coo, Z, XOUT = X_nodes, YOUT = Y_nodes, missing=!Values.F_NAN )



==== tiles ====


Sometimes it might be preferable to have just colored tiles around the nodes of your data. Here is an example how to do it:
The coordinate system must be defined - you may use the contour command as above with the /nodata keyword

  ; field sizes
  Nx = n_elements(x_coo)
  Ny = n_elements(y_coo)
  ; predefine neighbour nodes
  x2 = x_coo[0]
  x1 = x2
  ; intermediate border
  x12 = (x1+x2)/2.
  for ix = 0, Nx-1 do begin
    ; shift through nodes and borders ...
    x0 = x1
    x1 = x2
    x01 = x12
    ; new node and border
    x2 = x_coo[min([ix+1,Nx-1])]
    x12 = (x1+x2)/2.
    ; predefine ...
    y2 = y_coo[0]
    y1 = y2
    y12 = (y1+y2)/2.
    for iy = 0, Ny-1 do begin
      ; shift through ...
      y0 = y1
      y1 = y2
      y01 = y12
      ; new node ...
      y2 = y_coo[min([iy+1,Ny-1])]
      y12 = (y1+y2)/2.
      ; value
      zij = data[ix,iy]  
      ; color index
      if zij le min_data then i_col = 0 else $
      if zij le max_data then i_col = fix(float(zij-min_data)/(max_data-min_data)*(N_iso-1)) $
                         else i_col = N_iso-1
      ; color from color index
      color = iso_colors[i_col]
      ; fill the rectangle
      polyfill, [ x01, x01, x12, x12 ] , [ y01, y12, y12, y01 ], color = color
    endfor ; iy
  endfor ; ix 

==== arrows ==== 

IDL can draw arrows (or vectors). The basic command is arrow which takes the starting point and end point of a (or several) vectors. It can take coordinates for the data-, normalized- or device system. 

  arrow, x0,y0, x1,y1, [, color=...][, thick=...][ ,/data | /normal | /device ][, ...]

Allthough it accepts the /polar keyword it is ignored. If you have vectors in polar coordinates you will need to convert them with cv_coord to cartesian coordinates (=rectangular).\\
If you have a non isotropic coordinate system with different scaling and lengths of the x-axis and y-axis you need to correct for this:
  * we want to have in x and y the same device coordinates for an amplitude a in u or v-direction
  * to achieve this we mutliply the v component of the vector with a factor fac:
  * => dx_dev = a * !X.S[1]*!D.x_size = a * fac * !Y.S[1]*!D.y_size = dy_dev
  * if we want to have dx_dev = dy_dev we need
  * fac = !X.S[1]*!D.x_size / (!Y.S[1]*!D.y_size) 
Within factor the ratio !X.S[1]/!Y.S[1] corrects for the different scaling of the axes and the ratio !D.x_size/!D.y_size corrects for the different lengths of the axes on the paper.

Example for N vectors (u[i],v[i]) at N positions x[i], y[i]

  factor = !X.S[1] / !Y.S[1] * !D.x_size / !D.y_size
  arrow, x, y, x+u, y+v*factor, /data
  
Similar with function velovect:

  plot, x,y, /nodata
  factor = !X.S[1] / !Y.S[1] * !D.x_size / !D.y_size
  velovect, u, v*fac, x, y, /overplot

Note that you first have to define the coordinate system with e.g. the plot command such that !X and !Y are correctly defined.






==== histogram ====

IDL knows one and two dimensional histogram. 

=== 1-D histogram ===

A one dimensional histogram is e.g. generated with

  binsize = 0.1
  histo = histogram( data , binsize=binsize, locations = classes , /nan )

the classes variable will contain the **lower borders** of the bins with the first bin starting at min(data). \\
(histogramm is made via %%histo=lonarr((max-min)/binsize+1)& for i do ii=(data[i]-min)/binsize)& data[ii] += 1& ...%%

A //quick and dirty// plot can be generated with psym=10:

  plot, classes, histo, psym=10

but as classes contains the lower borders of bins it is shifted a half binsize to the left and the plot begins and ends with an open box. Slightly better is to shift the classes and add zero values at begin and end.

  N_classes = n_elements(classes)
  classes_x = [ classes[0]-bin , classes, classes[N_classes-1]+bin ] + bin/2
  histo_x = [ 0 , histo, 0 ]
  plot, classes_x , histo_x , psym=10, $
    xtitle = 'class (unit)', $
    ytitle = 'freq.'

or you write an own plot-histo procedure ... or you ask jan :-) ...


A histogram with geometric spacing, i.e. each bin is by a factor q_bin larger than the previous.\\
In this case data points lower or equal to zero are not allowed:
  q_bin = 2.0
  log_binsize = alog10(q_bin)
  histo = histogram( alog10(data[where(data gt 0)]), binsize=log_binsize, locations = log_bins, /nan )
  classes = 10^log_bins

Plotting shall use then logarithmic scaling for the x-axes:

    plot, classes , histo , psym=10, $
      xtitle = 'class (unit)', $
      ytitle = 'freq.', $
      /xlog


=== 2-D histogram ===

Assume you have two parallel data series data_1 and data_2 and want a joint histogram :

  ; set bin-size, min and max and vector with bin-borders for data set 1
  bin_size_1 = 0.1 ; give an adequate value
  min_data_1 = min(data_1, /nan )
  max_data_1 = max(data_1, /nan )
  N_bins_1 = fix((max_data_1-min_data_1)/bin_size_1)+1
  ; output of histo refers to lower border of bins => add 1/2binsize to center them in bin
  bins_1 = min_data_1 + (findgen(N_bins_1) + 0.5) * bin_size_1 
   
  ; set bin-size, min and max and vector with bin-borders for data set 2
  bin_size_2 = 100. ; give an adequate value
  min_data_2 = min(data_2, /nan )
  max_data_2 = max(data_2, /nan )
  N_bins_2 = fix((max_data_2-min_data_2)/bin_size_2)+1
  ; output of histo refers to lower border of bins => add 1/2binsize to center them in bin
  bins_2 = min_data_2 + (findgen(N_bins_2) + 0.5)*bin_size_2
   
  ; make 2d Histogramm
  histo_2d = hist_2d( $
    data_1[*] , data_2[*] , $ ; using [*] allows you to put together multidimensional data as eg. images ...
    bin1 = bin_size_1, min1 = min_data_1, max1 = max_data_1, $ ; if there are nans in data idl complains bin1 is not allowed => provide min and max to circumvent
    bin2 = bin_size_2, min2 = min_data_2, max2 = max_data_2 $
    )
   
  ; plot it
  contour, histo_2d, bins_1, bins_2, $
    Nlevels = 20, $ ; some value
    /follow, $ ; annotate iso lines
    ; /cell_fill, $ ; or use color shading
    xtitle = 'data 1', $
    xstyle = 1, $
    ytitle = 'data 2', $
    ystyle = 1



=== 2-D histogram vs time of day ===


Assume you have a time series and want a 2-D-histogram with values versus hour of day

  ; get hour from time (given as julian date)
  caldat, time, month, day, year, hour, minute, second
  hour_float = hour + (minute+second*60.)*60.
  ; alternatively
  ; hour_float = ((time-0.5-long(time-0.5)))*24.0
  
  ; get data 
  data = fltarr( N_time )
  ...
  
  ; define bin sizes
  bin_data = 100.0
  bin_hour = 1.0
   
  ; construct classes
  data_min = min( data )
  data_max = max( data )
  N_data_class = floor((data_max-data_min)/bin_data)+1
  class_data = data_min + findgen(N_data_class)*bin_data + bin_data/2
  
  ; make 2d Histogramm
  histo_2d = hist_2d( data , hour_float , bin1 = bin_data, bin2 = bin_hour )
  
  


===== 3-D grafics =====

==== 2 dim surface ====

To display a 2 dim function f(x,y) as 3-dim surface as a grid mesh use %%Surface%% %%shade_surf%%:

<code idl>
  ; generate function data
  Nx = 100
  Ny = 100
  ; x and y = -1..1
  x = 2*findgen(Nx)/Nx-1
  y = 2*findgen(Ny)/Ny-1
  ; N*N arrays with x and y - coordinates
  xx = x  #  replicate(1,Ny) ; # is matrix product, this gives a Nx*Ny matrix with x increasing along x direction
  yy = replicate(1,Nx)  #  y ; # is matrix product, this gives a Nx*Ny matrix with y increasing along y direction
  ; 2 dim radius array
  rxy = sqrt( xx^2 + yy^2 )
  ; a nice peaked function 
  fxy = cos(rxy*2*!pi*3) * exp(-(rxy/0.3)^2)

  ; plot in 3D ...
  surface, fxy,x,y, xtitle = 'x', ytitle='y', ztitle='z'

</code>
IF you want to have a solid, shaded surface use %%shade_surf%% instead:
<code idl>
  shade_surf, fxy,x,y, xtitle = 'x', ytitle='y', ztitle='z'
</code>


==== 3dim iso surface ====

To plot a 3dimensional isosurface of a 3D-field use %%shade_volume%% to get the surface and %%polyshade%% to draw it:
<code idl>
  ; calucalte isosurface from data = fltarr( Nx, Ny, Nz )
  shade_volume, data, iso_value, vertices, indices

  ; convert vertices given in index coordiantes to data coordinates
  vertices[0,*] = xcoo[ vertices[0,*] ]
  vertices[1,*] = ycoo[ vertices[1,*] ]
  vertices[2,*] = zcoo[ vertices[2,*] ]
  
  ; set up 3D transformation 
  SCALE3, XRANGE=[xmin,xmax], YRANGE=[ymin,ymax], ZRANGE=[zmin,zmax]

  ; draw x-y-z axes and ...
  ; ... use data from bottom of volume for basic zero plane
  contour, $
      data[ * , * , 0 ], $
      xcoo, $
      ycoo, $
      zvalue = 0.0, $ ; relative position on z-axis
      nlevels = 255, $
      /fill , $
      /cell_fill, $ ; is necessary otherwise contours are not closed or reach diagonally over the plane ...
      xtitle='x (km)', $
      ytitle='y (km)', $
      ztitle='z (km)', $
      zrange=[ zmin, zmax ], $ ; otherwise z will be rescaled to to the range of data[*,*,0]
      /t3d, $ ; use the 3dim transformation defined above
      charsize=2.5

  ; draw iso surface
  img = polyshade( vertices , indices , /T3D )
  
  ; save to png file
  save_png, 'iso_surf.png', img

</code>

If no device was established %%img%% will contain the image data. If the **z-buffer** is used (%%set_plot,'z'%%) variable img will contain a dummy value and the image can be get by use of **%%tvrd()%%**. If you are using the z-buffer several subsequent calls to %%polyshade%% and %%contour%% (with keyword %%/overplot%%) can be made.


==== 3 dim data ====

To display 3 dimensional data with semi transparent coloring use function %%voxel_proj%%.

There are several more functions see [[http://www.physics.nyu.edu/grierlab/idl_html_help/funclisting2.html#wp1184252|3D Visualization]] in the [[http://www.physics.nyu.edu/grierlab/idl_html_help/funclisting.html|Categorial list of functions]].


==== principles ====

IDL uses the system matrix %%!P.T%% to describe the 3D-projection. It can be manipulated with the functions 

  * %%surfr%% to define the viewing angles as rotation angles aroud the axes
  * %%scale3%% to define scaling along x,y,z and z-axis and to define the viewing angles
  * %%t3d%% reset, set and manipulate the transformation matrix and its parts

There is an article [[http://en.wikipedia.org/wiki/Perspective_projection|Perspective projection]] in the english wikipedia.
The basics are: 

We want to map an **object point** in object space **A**=(Ax,Ay,Az) to an **image point** on the screen **B**=(Bx,By,Bz) where Bz is its position relative to the screen (behind or in front - important for hidden lines).

The observer or camera is at C=(Cx,Cy,Cz), its viewing direction is defined by the 'eulerian' angles t=(tx,ty,tz) and is position relative to the screen by E=(ex,ey,ez) (i.e. E is given in the coordinate system of the camera). 

The following steps are performed:

  * %%D' = A-C%% position of A relative to the camera 
  * %%D = T*D'%%  rotate D into the camera system with rotation matrix T = f(tx,ty,tz) = [ [ Txx, Txy, Txz ] , [ Tyx , ... ] , ... ] 
  * %%B = (D-E)*ez/dz%%  project D onto the screen plane
  
Idl system matrix !P.T is a 4*4 Matrix describing the operation %% B' = !P.T * (A-C)%%

  Bx'    =  (  Txx  ,  Txy  ,  Txz    ,  -ex   )  *  (  Ax-Cx  )
  By'    =  (  Tyx  ,  Tyy  ,  Tyz    ,  -ey   )  *  (  Ay-Cy  )
  Bz'    =  (  Tzx  ,  Tzy  ,  Tzz    ,  -ez   )  *  (  Az-Cz  )
  wz'    =  (   0   ,   0   ,  -1/ez  ,   1    )  *  (    1    )

or something like that ... (there is something wrong with -ez and -1/ez)

The image coordinates are calculated then by Bx = Bx'/wz'.
  
The rows of rotation matrix T can be interpretad as the unit vectors describing the orientation of the screen in the object coordinate system.

===== 3D into pdf =====

PDF supports in principle 3D grafics - i.e. the PDF standard describes a format to incorporate 3D models into PDF. This format is supported by the ADOBE reader since version 8 or 9, by other readers usually not. Examples can be found in this [[http://compgroups.net/comp.lang.idl-pvwave/making-3d-pdf-from-idl/2089663|idl thread of compgroups.]]\\
[[http://dx.doi.org/10.1111/j.1365-2966.2011.19900.x]]\\
[[http://dx.doi.org/10.1038/nature07609]]\\
[[http://dx.doi.org/10.1126/science.1220574]]\\
A tutorial how to produce such 3D grafics can be found at [[http://www.astrobetter.com/tutorial-for-embedding-3d-interactive-graphics-into-pdf/| www.astrobetter.com]]. Three 3 steps are suggested:
  - produce a 3D representation of the data in [[http://de.wikipedia.org/wiki/Wavefront_OBJ|.obj code]] i.e. a (ascii) text representation
  - convert the .obj files into [[http://de.wikipedia.org/wiki/Universal_3D|universal 3D]] a XML standard for hierarchical 3D data
  - generate a PDF file with [[http://www.latex-project.org/|LaTeX]] and PDFLaTeX and embed the .u3d data by using the [[http://www.ctan.org/pkg/movie15|movie 15 package]]



===== Time/Date =====

current time as  string in the standard form of your system 
  t = SYSTIME()

current time in UTC as string ...
  t = SYSTIME( /UTC )

current time in [[https://en.wikipedia.org/wiki/Julian_day|Julian Day]] (Days (in fractions) since November 24, 4714 BC, 12:00:00 UTC)
  t = SYSTIME( /JULIAN, /UTC )
Actually /julian alone gives not UTC - you have to add the flag /UTC


current time in [[https://en.wikipedia.org/wiki/Unix_time|unix epoch]] (seconds since 1.1.1970 00:00:00 UTC - watch out for [[https://en.wikipedia.org/wiki/Year_2038_problem|Jan 19, 2038]])
  t = SYSTIME(/seconds)
this time contains fractions of seconds, can be used to find out how long program steps take:

  t0 = SYSTIME(/sec)
  ... do something ...
  t1 = SYSTIME(/sec)
  print, 'this took', t1-t0,'sec', f='(A,f0.3,A)'

current time as string formatted humand readable downto milliseconds:
  t = STRING( SYSTIME(/JULIAN), FORMAT='(C(CDI02,".",CMoI02,".",CYI04,"/",CHI02,":",CMI02,":",CSF06.3))')
Actually SYSTIME(/JULIAN) makes only steps of full seconds, accordingly the fractional seconds will remain zero.
But we can use the (unix-)seconds output, convert it to julian and print:

  sec_per_day = 24d*60d*60d
  ref_day = julday(1,1,1970,0,0,0)
  t = STRING( ref_day+SYSTIME(/SEC)/sec_per_day, FORMAT='(C(CDI02,".",CMoI02,".",CYI04,"/",CHI02,":",CMI02,":",CSF06.3))')



==== Julian day ====

Julian date = time in days since 1. Januar −4712 (4713 BC) 12:00 TDT ( = UTC + 1min 5.184sec )\\
see [[http://en.wikipedia.org/wiki/Julian_date]] or [[http://de.wikipedia.org/wiki/Julianisches_Datum]]

Julian day is a continous time variable and can easily be used to 
  * calulate time differences (dt = t2-t1 in days, dt=(t2-t1)*24. in hours ... )
  * print formated times (see below)
  * annotate grafics 

julian day is //the IDL format for time//.
It is stored as a double variable counting the days.

Be aware that zero fractional part refers to 12:00 (as the start time of julian day is 12:00)

=== create and decode time ===

To create a time in julian date use function julday:

  time = julday( month, day, year, hour, minute, second )
 
  * You can omit the time (hour, minute, second). 
  * If you provide at least the hour you get a double precision variable, if not, you get an integer variable of type long.
  * If one of the input variables are arrays you get an array. 
  * If these arrays have different size IDL uses the smallest size and gives no(!) warning. 
  * Zero input values for all except year are allowed
  * If year is zero you get an error message //"there is no year zero in the civil calendar"//
  * Negative input values are allowed.
  * You can provide strings as input, but be aware that empty strings are interpreted as zero.

To calculate a date use procedure caldat :

  caldat, time, month, day, year, hour, minute, second

  * if time is an array output month, day, etc will also be arrays



=== print formatted time ===

print formatted date from julian day in //two steps//:
  CALDAT, t, M,D,Y,HH,MM,SS
  PRINT, D,'.',M,'.',Y,' ',HH,':',MM,':',SS,FORMAT="(i2,a,i2,a,i4,a,i2,a,i2,a,i2)"
or more elegant in //one step// for DD.MM.YYYY hh:mm:ss :
  PRINT, t, FORMAT='(C(CDI02,".",CMOI02,".",CYI04," ",CHI02,":",CMI02,":",CSI02))'

or for YYMMDDhhmmss :
  PRINT, t, FORMAT='(C(CYI02,CMOI02,CDI02,CHI02,CMI02,CSI02))'

you can also have fraction of a second for hh:mm:ss.ttt i.e. seconds downto 1/1000.sec use:
  PRINT, t, FORMAT='(C(CHI02,":",CMI02,":",CSF06.3))'



the function STRING can be used as well; for DD.MM.YYYY hh:mm:ss use:
  s = STRING( t, FORMAT='(C(CDI02,".",CMOI02,".",CYI04," ",CHI02,":",CMI02,":",CSI02))' )
for YYMMDDhhmmss use:
  s = STRING( t, FORMAT='(C(CYI02,CMOI02,CDI02,CHI02,CMI02,CSI02))' )

The format specifiers are construced as:
  C<field identifier><format specifier>[<leading zero>]<number of digits>
  * <field identifier> can be D for day, MO for month, Y for year, H for hour ...
  * <format specifier> can be F for float, I for integer, A for character
  * <leading zero> give a zero for leading zeros
  * <number of digits> 2 makes a field with a width of 2 digits/characters

(see Help for [[http://www.physics.nyu.edu/grierlab/idl_html_help/files24.html#wp2823814|"format codes" -> "C() format codes"]] )



=== Rounding ===

The string-formatting routine does //not round// the fields, but instead delivers the floor of each field.
This may result in e.g. <hh-1>:59 if you ought to have hh:00 but have a 1/10000s missing.
To avoid this add half of the smallest unit you do print and round times:
  * for hh:mm    use %%print, round(time*minperday+0.5)/minperday, format='(C(CHI02,":",CMI02))'%%
  * for hh:mm:ss use %%print, round(time*secperday+0.5)/secperday, format='(C(CHI02,":",CMI02,":",CSI02))'%%
  * for hh:mm:ss.ttt use  %%print, time, format='(C(CHI02,":",CMI02,":",CSF6.3))'%% (fractions of a second are rounded ???)
with 
  minperday = 24d*60d = 1440.
  secperday = 24d*60d*60d = 86400.


=== lowest resolvable unit ===

Julian day comes in [[https://en.wikipedia.org/wiki/Double-precision_floating-point_format|type DOUBLE-Precision]] which has a mantissa of 52 bits. The smallest resolvable step is thus dt_rel=1/2^52=2.2E-16. As we have currently (2019) Julian day number t=2458695 this is dt_min=dt_rel*t*sec_per_day=4.7E-5sec=0.047msec (sec_per_day=24.*60*60=86400). You can test this with this code (:-/):

  spd=86400d&t0=julday(1,1,1970,0,0,0)&twt=1e-6&for i=0,99 do begin&t1=systime(/sec)/spd+t0&wait,twt&t2=systime(/sec)/spd+t0&print,(t2-t1)*spd&end

It will most of the time print 0.0000 and sometimes 4E-5. If you lower twt (t_wait) it will always print 0.0000



=== convert time formats ===


**'human' date to julian date**
  t = JULDAY(Month, Day, Year [,hour[,minute[,second]]] )
if you omit hour, minute, second you will get type long, otherwise double

**string to julian date** 
for a string of the form %%DD.MM.YYYY hh:mm:ss:%%
  ; DD.MM.YYYY hh:mm:ss
  ; 0123456789 12345678
  ;   JULDAY(   month      ,   day        ,   year       ,    hour       ,    minute     ,    second      ) 
  t = JULDAY( strmid(s,3,2), strmid(s,0,2), strmid(s,6,4), strmid(s,11,2), strmid(s,14,2), strmid(s,17,2) )



**Time differences** can be easily calculated - but be aware that julian date starts at 12:00 - if you want **display** time differences as time strings you must subtract a half day:
  PRINT, t1-t2 - 0.5 , FORMAT='(C(CHI02,":",CMI02,":",CSI02))'

 


To **annotate a plot axis** with formatted time strings just set an appropriate format 
specifier in the plot command:
  plot, time, data, ..., xtickformat = '(C(CHI02,":",CMI02))', ...
where time is given in julian day ...

==== Unix time (seconds since ...) ====

Unix time are the seconds since 1.1.1970 00:00:00 UTC \\
Sometimes it is also called //Epoch// (see [[http://en.wikipedia.org/wiki/Unix_time]])

Originally t_unx is defined as long integer - this leads to the [[http://en.wikipedia.org/wiki/Year_2038_problem|year 2038 Problem]] 


Julian day (t_jul) to unix time (t_unx):
  secperday = 24*60*60d
  unx_start = JULDAY( 1, 1, 1970, 0, 0, 0 ) ; this is equal to 2440587.5
  t_unx = (t_jul - unx_start) * secperday

and vice versa:
  t_jul = t_unx/secperday + unx_start

==== Other time formats ====

If you have a time t_xxx as %%tics%% since %%t_ref%% (e.g. (fractional)number of 7min intervals since 3.4.1722 13:21 (just to show that it can be really weird :-) ) then
define:
  tics_per_day = min_per_tic / min_per_day  or  sec_per_tic/sec_per_day  etc...
  t_ref = JULDAY( mo_ref, dd_ref, yy_ref,  HH_ref, MM_ref, SS_ref )

convert julian day (t_jul) time into this exotic format (t_xxx)
  t_xxx = (t_jul - t_ref) * tics_per_day
and vice versa:
  t_jul = t_xxx/tics_per_day + t_ref


Examples are:
  * **excel time 1** = days since 1.1.1900 (nineteen hundred) \\ (excel has two internal timeformats, you can switch somewhere inside excel between both)
    * t_ref = julday( 1,1,1900, 0,0,0 )
    * tics_per_day = 1.0
    * typical values (june 2023) are in the order of 45087 (45tsnd) days
    * Excel cannot deal with negative times, i.e. dates before t_ref are not valid
    * if you export excel tables as csv you might get this as time variable
  * **excel time 2** = days since 1.1.1904 (nineteen O four! )
    * t_ref = julday( 1,1,1904, 0,0,0 )
    * tics_per_day = 1.0
    * typical values (june 2023) are in the order of 43627 (43tsnd) days
    * if you export excel tables as csv you might get this as time variable
    * you can switch somwhere inside excel between both time systems - be careful
    * Excel cannot deal with negative times, i.e. dates before t_ref are not valid
    * Excel has a known bug in the year 1900 it thinks it errornously assumes it was a leap year => all excel dates before 1.3.1900 might be wrong (weekdays, differences etc.)
  * **unix time** = seconds since 1.1.1970
    * t_ref = julday( 1,1,1900, 0,0,0 )
    * tics_per_day = 86400   (=24*60*60)
    * typical values (june 2023) are in the order of 1686589673 = 1.67E9 sec
    * this ist the internal time format on unix systems and derivatives (e.g. android)
    * unix time may experience an overflow or a roll over in 2038 when range of long32 is exceeded. This may lead to problems in the computer world.
  * **days since 1.1.2000**  ( time format in RPG data ? )
    * t_ref = julday(1,1,2000,0,0,0)
    * tics_per_day = 1.
    * typical values (june 2023) are in the order of 8563 days
  * **seconds since 1.1.2000**  ( time format in Juelich radiation data )
    * t_ref = julday(1,1,2000,0,0,0)
    * tics_per_day = 86400.
    * typical values (june 2023) are in the order of 739905199 = 7.3990510e+008 = 740Mio seconds
  


==== day of year, hour of day ====

The day of the year (%%doy%%) - sometimes incorrectly named 'Julian day' - can be calulated from a Julian date by:
<code idl>
CALDAT, t_julian, month, day, year
t_doy = floor( t_julian - JULDAY( 1,1,year ) ) + 1
</code>
Note: Day of year starts with 1 on January first => add one.
The simplified formula %%doy = t_julian - floor(t_julian/days_per_year)*days_per_year%% - with %%days_per_year=365.25%% does not work - why ?

The hour of the day (%%hod%%) is the fractional part of the date as julian day shifted by half a day/twelve hours:
<code idl>
t_hod = ( t_julian-0.5 - floor( t_julian-0.5 ) ) * 24.
</code>

or alternatively and more intuitive, but slower
<code idl>
CALDAT, t_julian, month, day, year, hh, mm, ss
t_hod = hh+( mm + ss/60.)/60.
</code>




===== interpolation =====

IDL has some predefined functions for interpolation:
  * [[https://www.nv5geospatialsoftware.com/docs/INTERPOL.html|interpol]] interpolates a 1dim vector to new abscissa values - either on a regular grid or with given abscissas
  * [[https://www.nv5geospatialsoftware.com/docs/interpolate.html|interpolate]] can do bilinear interpolation on two-dimensional regular or irregular grids
  * [[https://www.nv5geospatialsoftware.com/docs/bilinear.html|bilinear]] can do bilinear interpolation on regular grids
  * [[https://www.nv5geospatialsoftware.com/docs/krig2d.html|krig2d]] does 2 dimensional Kriging i.e. a distant weighted interpolation from neighboring points.

The interpolate routine takes only floating point index values for the interpolation. If you want to interpolate to certain nodes x_ipol and y_ipol you have to calculate the respective indices. This can be done by the interpol routine as follows. Assume we have a data array z_org = fltarr(Nx_org,Ny_org) with coordinates x_org=fltarr(Nx_org) and y_org=fltarr(Ny_org) and want interpolate it to z_ipol = fltarr(Nx_ipol,Ny_ipol) with coordinates x_ipol=fltarr(Nx_ipol) and y_ipol=fltarr(Ny_ipol):

  ; bilinear interpolation 
  ix_ipol = interpol( findgen(Nx_org) , x_org , x_ipol )
  iy_ipol = interpol( findgen(Ny_org) , y_org , y_ipol )
  z_ipol = interpolate( z_org , ix_ipol , iy_ipol , /grid  )

If you want to interpolate only over certain dimensions they must be the last ones. You may use transpose to rearrange. Lets assume you have an array depending on x, y and time, i.e. z=fltarr(Nx,Ny,Nt) and you want to interpolate to certain x and y coordinates you can do this as follows:

  z_ipol = transpose( interpolate( transpose(z_org,[2,0,1]) , ix_ipol , iy_ipol , /grid  ), [1,2,0] )



===== geographical coordinates =====

IDL can transform geographical coordinates into cartesian  x,y coordinates and vice versa. 
To do that you have to define your coordinatesystem with %%MAP_PROJ_INIT%%. To calculate your projected coordinates %%(X,Y)%% from %%longitude%% and %%latitude%% use %%MAP_PROJ_FORWARD%% and to get %%longitude%% and %%latitude%% for a %%(X,Y)%% pair use %%MAP_PROJ_INVERSE%%.

Example for the universal transverse mercator projection (UTM) for West Germany (UTM zone 32 with central meridian at 9°East):

Init the transformation:
<code idl>
  map_struct = MAP_PROJ_INIT( $
    101 , $ ; UTM
    CENTER_LONGITUDE = longitude , $
    )
</code>
The result %%map_struct%% contains information of the projection and must be handed to the functions for the caluclation of coordinates. The parameter %%CENTER_LONGITUDE%% must be given to define which UTM grid to be used (read [[http://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system|wikipedia]] if you want to understand the details). It should be a longitude close to the locations you are interested in. Following the UTM definition IDL will round it to a multiple of 6° relative to 3° (i.e. %%CENTER_LONGITUDE=3+n*6%% with n an integer).
Calcualted UTM coordinates will be conform to the UTM specifaications but the direction of UTM-northing may deviate by several degrees from true north. If you need a coordinate system directed exactly to north  at your location use instead:

<code idl>
  map_struct = MAP_PROJ_INIT( $
    109 , $ ; transverse mercator
    DATUM = 8 , $ ; WGS84
    CENTER_LONGITUDE = longitude , $
    FALSE_EASTING = 500000. $
    )
</code>

Calculate some UTM values from longitude %%long%% and latitude %%latt%%:
<code idl>
  UTM = MAP_PROJ_FORWARD( long , latt , map_structure=map_struct  )
</code>
where long and latt can be arrays (fltarr(N)) and UTM will be of type fltarr(2,N) containing the calculated UTM coordinates.

The inverse: calculate some longitude and latitude values from UTM=fltarr(2,N)
<code idl>
  LonLat = MAP_PROJ_inverse( UTM, map_structure=map_struct )
</code>




A list of the ID's of the projections can be found at [[http://www.harrisgeospatial.com/docs/map_proj_init.html#M_824365677_981059|MAP_PROJ_INIT]] help article.


You can also generate a map and draw into this map by using ''MAP_SET'' :

<code idl>

  center = [ center_lat , center_lon ]
  d_ll = 30.
  map_box = [ center - [1,1]*d_ll , center + [1,1]*d_ll ]

  MAP_SET, /Azimuthal, /isotropic, $
    center_lat , center_lon , 0 , $
    limit = [ lat_a, lon_a, lat_b, lon_b ]
   
  ; project an image to the map with img=bytarr(Nx,Ny)
  read_jpeg, '/data/sfb1211/A01/srtm/other/images/Envisat_global_land_cover_map.jpg', img
  latmin_jpg = -94. &  latmax_jpg = +106. & lonmin_jpg = -180. & lonmax_jpg = +180.
  ; position to put image
  start_x = lonmin_map
  start_y = latmin_map
  ; wrap image - MAP_image can deal only with 2dim array - we need to do it channel by channel
  img_mapped_red = MAP_IMAGE( reform(img[0,*,*]) , startX, startY, lonmin=lonmin_jpg, lonmax=lonmax_jpg, latmin=latmin_jpg, latmax=latmax_jpg )
  img_mapped_grn = MAP_IMAGE( reform(img[1,*,*]) , startX, startY, lonmin=lonmin_jpg, lonmax=lonmax_jpg, latmin=latmin_jpg, latmax=latmax_jpg )
  img_mapped_blu = MAP_IMAGE( reform(img[2,*,*]) , startX, startY, lonmin=lonmin_jpg, lonmax=lonmax_jpg, latmin=latmin_jpg, latmax=latmax_jpg )
  ; put channels togehter to one image
  img_mapped_rgb = transpose( [ [[img_mapped_red]], [[img_mapped_grn]], [[img_mapped_blu]]], [2,0,1] )
  ; put image on canvas
  TV, img_mapped_rgb, startX, startY, true=1


  ; draw coastlines, continents and Country borders in high resolution 
  MAP_CONTINENTS, /COASTS, /CONTINENTS, /COUNTRIES, /HIRES

  ; rivers - idl draw allways all rivers - an hierarchy as in GMT is not there (https://docs.generic-mapping-tools.org/6.3/pscoast.html#i)
  MAP_CONTINENTS, /RIVERS, /HIRES, color = 'ff0000'x
  
  ; draw a grid, make nice black and white axes around it
  MAP_GRID , /box_axes
  
  ; plot and annotate some data (lat=fltarr(N), lon=fltarr(N), name=strarr(N))
  oplot, lon, lat, psym=4
  for i = 0, n_elements(lat)-1 do xyouts, lon[i], lat[i], name[i]
  
</code>
===== wind directions =====

IDL's **atan function** can be fed with x and y component of a vector to get the direction angle from -pi to +pi = -180° to +180°. 

  phi = atan( x, y )

It performs the complete angle determination based on the signs of x and y, i.e. something like

  function atan_2, y, x
    ; determination of azimuth angle from x and y components
    ; this is principally done if you call atan( y,x )
  if  x eq 0. then begin
    if y gt 0. then phi =  !pi/2. else $
    if y lt 0. then phi = -!pi/2. $
               else phi =  0.0
  endif $; x eq 0
  else begin ; x ne 0
    phi0 = atan( y/x )
    if (y gt 0. ) and (x gt 0.) then phi = phi0       else $ ; first  quadrant
    if (y le 0. ) and (x gt 0.) then phi = phi0       else $ ; second quadrant
    if (y lt 0. ) and (x lt 0.) then phi = phi0 - !pi else $ ; third  quadrant
    if (y gt 0. ) and (x lt 0.) then phi = phi0 + !pi        ; fourth quadrant
  endelse 
  return, phi
  end ; function 

As the //wind directon azimuth angle// is differently defined than mathematical azimuth angle i give here the ultimative conversion laws between **wind direction** %%azi%%, wind speed %%spd%% and **vector components** %%u%% (W-E) and %%v%% (S-N):

  ; spd and azi  to  u and v
  u = - spd * sin(azi*!dtor)
  v = - spd * cos(azi*!dtor)

  ; u and v  to  spd and azi
  azi = atan(  -u , -v ) / !dtor
  if azi lt 0 then azi += 360.
  spd = sqrt( u^2 + v^2 )

Factor ''!dtor'' converts **d**egrees **t**o **r**adians, it is of type FLOAT (=single pecision) for DOUBLE precisions use ''!CONST.DtoR''

Examples: 
  * North wind: azi = 0° => u = 0 (no W-E component), v = -spd (wind is blowing from north to south)
  * East wind: azi = 90° => u = -spd (wind is blowing from East to West), v = 0 (no N-S component)
  * West wind: azi =270° => u = +spd (wind is blowing from West to east), v = 0 (no N-S component) 

==== Wind direction error ====


The error can be estimated with gaussian error propagation:

<m>
  \phi = arctan(-u/-v)
       = arctan(x)  \ ; \  
  x = {-u}/{-v}
</m>

<m>
  \Delta \phi^2 = ({d\phi}/du \Delta u)^2 + ({d\phi}/dv \Delta v)^2
               = ({d\phi}/dx dx/du \Delta u)^2 + ({d\phi}/dx dx/dv \Delta v)^2
</m>     
          
if we consider <m>{d\phi}/dx=1/{1+x^2}</m> we get 

<m>
  \Delta \phi^2 = 1/(u^2+v^2)^2 ( (v \Delta u)^2 + (u \Delta v)^2 )
</m>


===== percentiles =====

(or quantiles)

IDL can give you the median but not any other arbitrary percentile. So do it yourself:

  ; define 5 percentiles
  N_perc = 5
  p1 = 0.05 ; 5% and 95%
  p2 = 0.25 ; 25% and 75% = lower and upper quartiles
  ; add the median and mirror p1 and p2:
  p_fac = [ p1 , p2 , 0.5, 1-p2 , 1-p1 ]
  
  ; determine the percentiles
  percentiles = (data[ sort[data] ])[ p_fac*(n_elements(data)-1) ]

This way +/-infinity and Nan values are included. IDL sorts them in the following way: 
  data[sort[data]] = [[-inf.values], [finite.values], [+inf.values], [Nan.values] ]
You have to count the respective non-finite values and include it in the index calculation:
  N_neg_inf = long( total(finite(data,/infinity,sign=-1)) )
  N_pos_inf = long( total(finite(data,/infinity,sign=+1)) )
  N_nan     = long( total(finite(data,/nan             )) )
  percentiles = (data[ sort[data] ])[ N_neg_inf + p_fac*(n_elements(data)-1-(N_neg_inf+N_pos_inf+N_nan)) ]
An implementation as function can be found in /home/hatpro/idl/lib/f_percentile.pro

There is no exact method to determine percentiles. 
Imagine a dataset with 50% of the values lower than a vlaue y1 and the other 50% larger than a value y2 with a large gap between y1 and y2. 
  * Where ist the median in this case? Is it y1, y2 or somewhere in between?
  * Is interpolation between y1 and y2 allowed?
  * In the same way any percentile may not be exactly determined
  * [[http://en.wikipedia.org/wiki/Percentile|wikipedia]] gives a method for interpolation between these values...

===== random numbers =====

Random numbers can be generated with 

  x = RANDOMU( seed )
  x = RANDOMU( seed, N )
for a single number or an array with N elements of **U**niform distrinbuted random numbers in  [0,1)\\
Expectation value is thus μ=0.5 and standard deviation σ=3.46.\\
(Uniform disributed radnom numbers have an expectation value μ=(a+b)/2 and standard deviation σ=sqrt(12)*(b-a) if they lie in the range [a,b))

  x = RANDOMN( seed )
  x = RANDOMN( seed, N )
for as single number or an array with N elements of **N**ormal distrinbuted random numbers with expectation value μ=0 and standard deviation σ=1.0. There is theoretically no upper and lower border for these random numbers but |x-µ| will probably never be much larger than 4σ.

**BUT** there is no //Zufall// in computers !\\
All these routines produce //pseudo random numbers//\\
The random number is //calculated// from seed !\\
With every random number a new seed is returned for the next random number ...\\
The following code will allways produce the same sequence of ten //"radnom"// numbers:

  seed = 0
  for i =0, 10 do PRINT, RANDOMU(seed)

(you may also omit %%seed = 0%% as idl sets %%seed%% than internally to 0) 

To generate some randomness you can use the time:
  seed = SYSTIME(/seconds)
  for i =0, 10 do PRINT, RANDOMU(seed)

See [[http://www.nrbook.com/a/bookcpdf.php|Numerical recipies]] chapter 7 or [[http://de.wikipedia.org/wiki/Zufallsgenerator|wikipedia]]. 


===== binary-, octal-, decimal-, hexadecimal - numbers  =====
==== printing ====

every number can be displayed in many different ways. In IDL this canbe done with formatted output. \\
Let e.g. x = 16(decimal):
  * binary:  print, x, format='(B)' -> 10000
  * octal: print, x, format='(O)' -> 20 
  * decimal: print, x, format='(I)' -> 16
  * hexadecimal: print, x, format='(Z)' -> 10
  * roman: print, ääääää ... ok that does not work - but it should look like XVI

The principle of all these number representations (except for the roman) is [[http://en.wikipedia.org/wiki/Positional_notation|positional notation (wikipedia)]]:\\
A number written as a sequence of digits A<sub>i</sub>:\\
x = A<sub>3</sub>A<sub>2</sub>A<sub>1</sub>A<sub>0</sub>\\
is mathemtically the sum:\\
x = A<sub>3</sub>*P<sup>3</sup> + A<sub>2</sub>*P<sup>2</sup> + A<sub>1</sub>*P<sup>1</sup> + A<sub>0</sub>*P<sup>0</sup>\\
  * Where P is the base of the used number representation.
    * in the decimal system P=10
    * in the binary system P=2
    * in the octal system P=8
    * in the hexadecimalsystem P=16
    * ...
  * Digits A<sub>i</sub> can lie in the interval 0 to P-1.
    * in the decimal system (P=10) valid digits are 0..9
    * in the binary system (P=2) valid digits are 0..1
    * in the octal system (P=8) valid digits are 0..7
    * in the hexadecimalsystem (P=16) valid digits have the decimal values 0..15 - but since digits should be only one character wide the values 10 to 15 are replaced by the first six letters of the alphabet: 10->A, 11->B, 12->C, 13->D, 14->E, 15->F
  * Every position i in the written number is related to the base value P<sup>i</sup> 
    * the rightmost position in a decimal number has the base value 1 = 10<sup>0</sup>, \\ the position to the left has the base value 10 = 10<sup>1</sup>\\ the position to the left has the base value 100 = 10<sup>2</sup>\\ ...
    * the rightmost position in a binary number has the base value 1 = 2<sup>0</sup>, \\ the position to the left has the base value 2 = 2<sup>1</sup> \\ the position to the left has the base value 4 = 2<sup>2</sup>\\ ...
    * the rightmost position in a octal number has the base value 1 = 8<sup>0</sup>, \\ the position to the left has the base value 8 = 8<sup>1</sup> \\ the position to the left has the base value 64 = 8<sup>2</sup>\\ ...

Since computers do internally work with binary numbers it migth be sometimes useful to write binary numbers - but this creates easily very long numbers. That is the reason to use octal or hexadecimal numbers.  Their base is a power of 2 and thus compress some binary digits into one digit. In the octal system 3 binary postions are compressed into one (P<sub>oct</sub>=8=2<sup>3</sup>) and in the hexadecimal system 4 binary postions are compressed into one (P<sub>hex</sub>=16=2<sup>4</sup>):
  bin   oct      bin  oct ...   bin  hex(dec)    bin   hex(dec)
  000    0       100   4       1000   8  (8)    1100    C  (12)
  001    1       101   5       1001   9  (9)    1101    D  (13)
  010    2       110   6       1010   A (10)    1110    E  (14)
  011    3       111   7       1011   B (11)    1111    F  (15)


To convert a number x into any number system on the base P 'manually' (by mental arithmetic ?) divide the number by P and note the divisons rest (a number between 0 and P-1) as the rightmost digit. Continue the same with the result of the division and write the divisons rest to the left of your previous digit until you end up with zero ...


==== constants ====

To set IDL variables to certain binary, octal or hexadecimal values write the number in the desired representation as string followed by a letter inidcating the used number system:
  x = '10000'b
  x = '20'O
  x = '10'x

==== converting ====

To convert strings given in a number system to integer values use %%reads%% with an apropriate format specifier:
  ; binary
  s = '10000'
  x = 0L
  reads, s , x , format='(B)'
  print, x

  ; octal
  s = '20'
  x = 0L
  reads, s , x , format='(O)'
  print, x

  ; hexadecimal
  s = '10'
  x = 0L
  reads, s , x , format='(Z)'
  print, x


==== binary flags ====

In many cases the status of a system is delivered as flags in a binary number. Every position (i.e. bit OR flag) in the binary number is related to a certain status of a subsystem. A 0 represents //false// or //off// or //not active// etc. A 1 (or a set bit) represents vice versa a //true// or //on// or //active// etc. To check these //flags// it is convenient to use '//binary//' operators (here are operators meant dealing with binary numbers - the same word is used for operators between two numbers...):
  * c = a AND b : if bit at position i is set in a and b it will be set in c too and otherwise not:\\ a=1010<sub>bin</sub> b=0001<sub>bin</sub> => c=0000<sub>bin</sub>=0<sub>dec</sub>\\ a=1010<sub>bin</sub> b=0010<sub>bin</sub> => c=0010<sub>bin</sub>=2<sub>dec</sub>\\ ...
  * c = a OR b : if bit at position i is set in a or b it will be set in c and otherwise not:\\ a=1010<sub>bin</sub> b=0001<sub>bin</sub> => c=1011<sub>bin</sub>=11<sub>dec</sub>\\ a=1010<sub>bin</sub> b=0010<sub>bin</sub> => c=1010<sub>bin</sub>=10<sub>dec</sub>\\ ...
  * c = a XOR b : if bit at position i is set in a only or b only it will be set in c and otherwise not (exclusive or):\\ a=1010<sub>bin</sub> b=0001<sub>bin</sub> => c=1011<sub>bin</sub>=11<sub>dez</sub>\\ a=1010<sub>bin</sub> b=0010<sub>bin</sub> => c=1000<sub>bin</sub>=8<sub>dec</sub>\\ ...
  * c = NOT a : if bit at position i is not set in a it will be set in c and otherwise not:\\ a=1010<sub>bin</sub> => c=0101<sub>bin</sub>=5<sub>dez</sub>\\ ...

Dont mix these operators with //logical operators// 

To check whether in a status variable %%stat%% several flags are set define a %%mask%% and check by comparing %%stat%% with %%mask%%
using of the AND operator:
  ; define a mask to check for bits 2, 6, and 7
  mask = '00000100'b OR '01000000'b OR '10000000'b
  ; check
  IF (stat AND mask) ne 0 THEN BEGIN
    print, "some of the flags are set ..."
    ...
  ENDIF


==== hexadecimal representation ====

in rare cases it might be necessary to see the representation of a variable in memory. This can be done in the following way: 
  * copy content of variable into a long or long64 integer variable
  * print it formatted in hex, octal or binary
see [[https://www.nv5geospatialsoftware.com/Support/Self-Help-Tools/Help-Articles/Help-Articles-Detail/ArtMID/10220/ArticleID/19447/2723]]

For single precision variables we need a 4Byte=32Bit integer variable i.e. type long.
We print in Hexadecimal (formatcode Z), Octal (O) and Binary (B)
  x = 1.234
  xl = long(x,0)
  print, x, xl, xl, xl, f='(f," ",Z," ",O," ",B)'

For double precision variables we need an 8byte = 8*8bit = 64bit integer variable, i.e. Long64:
  x = 1.234d
  xl = long64(x,0)
  print, x, xl, xl, xl, f='(f," ",Z," ",O," ",B)'
  
This way we can see how NAN's and infinity are represented:
   print, long( !values.f_nan,0),      long64( !values.d_nan,0),      f='(" fnan=",Z,"  dnan=",Z)'
   print, long(-!values.f_nan,0),      long64(-!values.d_nan,0),      f='("-fnan=",Z," -dnan=",Z)'
   print, long( !values.f_infinity,0), long64( !values.d_infinity,0), f='("+finf=",Z," +dinf=",Z)'
   print, long(-!values.f_infinity,0), long64(-!values.d_infinity,0), f='("-finf=",Z," -dinf=",Z)'
will give:
   fnan=    7FC00000  dnan=      7FF8000000000000
  -fnan=    FFC00000 -dnan=      FFF8000000000000
  +finf=    7F800000 +dinf=      7FF0000000000000
  -finf=    FF800000 -dinf=      FFF0000000000000

The -NaN values is a bit strange as NaN is Not-a-Number but IDL also produces NaNs where the sign flag is set. 
===== Fourier Transform =====

A fourier transform is done with function FFT.
to plot the spectrum of a time series y running over a time interval of length T do the following:

  FT_y = FFT( y )
  
  N = n_elements( y )
  N2 = N/2
  i2 = 1+findgen(N2-1)
  freq = i2/N * 1/T
  FTy_cos = real_part( FT_y[i2] )
  FTy_sin = imaginary( FT_y[i2] )
  FTy_pow = FT_y[i2] * conj(FT_y[i2])
  
  plot, freq, FTy_cos, color='0000ff'x
  plot, freq, FTy_sin, color='ff0000'x
  plot, freq, FTy_pow

The power spectrum eventually needs some adjustment - depending on how you define it.

You can use FFT to calculate the auto-covariance and cross-covariance of two time series x and y:

  ft_x = FFT( x )
  ft_y = FFT( y )
  auto_cov_x = FFT( ft_x * CONJ(ft_x) , /INVERSE )
  cross_cov_xy = FFT( ft_x * CONJ(ft_y) , /INVERSE )

For long time series this is much faster than a classical sum( xi'*yi' ).

And you can do this also for two dimensional data sets x and y ... 


===== wavelet analysis =====

Wavelet analysis (WV) splits similar to [fast]fouriertransform (FFT) a dataseries into components at different scales (or wavelengths or frequencies). In difference to the FFT it gives also a localization of analysed scales, i.e. it gives information when or where in a dataseries a certain frequency had a a large amlitude.

Idl has a **Wavelet Toolkit** consisting of a number of graphical user interfaces and functions. It is described in the IDL help under [[http://www.physics.nyu.edu/grierlab/idl_html_help/introa2.html#wp996871|Introduction to the IDL Wavelet Toolkit]]. 

But you need to have a **license** for the toolkit to use the graphical user interface. We do not have the license! 

Nevertheless the basic functions are accessible without the license. In the help there is a [[http://www.physics.nyu.edu/grierlab/idl_html_help/ref2.html | list of the functions ]]. Documentation is rather poor sometimes it is helpful to study the source code of the functions you find in the //IDL-lib//%%/wavelet/source/%% directory (that is in our implementation in %%/opt/itt/idl/lib/wavelet/source/%%). It seems as if all the functions were written by Chris Torrence in 2000. He refers in his documentation to 

//[[http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%281998%29079%3C0061%3AAPGTWA%3E2.0.CO%3B2|Torrence, C. and G.P. Compo, 1998: A practical guide to wavelet analysis. Bull. Amer. Meteor. Soc., v. 79, pp. 61-78.]] [[http://dx.doi.org/10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2|DOI:10.1175/1520-0477(1998)079<0061:APGTWA>2.0.CO;2]]//



==== Wavelet transforms ====

Basically you can do two kinds of wavelet transformations, \\
the so called **discrete wavelet transform (DWT)** and the **continous wavelet transform (CWT)**: 
  * **DWT** is based on orthognal wavelets like the Haar or the Daubechie Wavelets. It splits the scale/space domain of the wavelet transform in separate regions with no overlapping information. It will deliver wavelet coefficients in the so called pyramidical order: two coefficents for the lowest scale, four for the second, ... and N coefficent for the N'th scale.
    * call with\\ %%wv = WV_DWT( data, Scaling, Wavelet, Ioff, Joff [, ... ] )%%
      * %%data = fltarray(N_data)%% input dataseries
      * %%scaling = fltarr(N_coeff)%% scaling coefficients of the wavelet, sometimes called the 'father wavelet', see //wavelet functions// below
      * %%wavelet = fltarr(N_coeff)%% wavelet coefficients, sometimes called the 'mother wavelet', see //wavelet functions// below 
      * %%Ioff = scalar integer%%, offset of the supprot for the scaling function, set to -N/2+2 to center it over each point in data
      * %%Joff = scalar integer%%, offset of the supprot for the wavelet function, set to -N/2+2 to center it over each point in data
      * [[http://www.physics.nyu.edu/grierlab/idl_html_help/ref8.html#wp1010303| further details on WV_DWT]]
  * **CWT** can be based on virtually any function (of course there are some restrictions) the best known are the morlet wave package or the mexican hat. The information received by the transformation has some overlap in the scale/space domain. It will deliver the wavelet coefficients as fltarr(N_data,N_scale)
    * call with \\ %%wv = WV_CWT( data, Family, Order [, ... ] )%%
      * %%data = fltarray(N_data)%%, input dataseries
      * %%Family = string%%, name of the wavelet function. Is expanded to %%WV_FN_%%//family//%%.PRO%%, see //wavelet functions// below
      * %%Order = float%%, parameter forwarded to the wavelet function, e.g. the number of zero intercepts of the morlet wavelet.
      * [[http://www.physics.nyu.edu/grierlab/idl_html_help/ref6.html#wp1009330| furhter details on WV_CWT]]
==== wavelet functions ====
 
IDL provides seven wavelet functions in directory %%/opt/itt/idl/lib/wavelet/source/%% with names
  wv_fn_coiflet.pro
  wv_fn_daubechies.pro
  wv_fn_gaussian.pro
  wv_fn_haar.pro
  wv_fn_morlet.pro
  wv_fn_paul.pro
  wv_fn_symlet.pro
These functions return a structure with information about the wavelet. The values of the wavelet funcion itself come in parameter wavelet.
Despite their similar names there are two different types of wavelet functions in IDL: //DISCRETE// and //CONTINOUS//:
  * **discrete** are the  coiflet, daubechies, haar and symlet wavelet
    * they can be used in the discrete wavelet transform %%WV_DWT%% only !
    * they are called by \\ %%wv_info = WV_FN_<name>( [Order, Scaling, Wavelet[, Ioff, Joff]] )%%
    * the 'mother wavelet' in array wavelet is returned in the third parameter without keyname
    * Scaling is //output// parameter containing the 'father wavelet' as array
    * wavelet is the wavelet function 
  * **continous** are the gaussian, morlet and paul wavelet
    * they can be used in the continous wavelet transform %%WV_CWT%% only !
    * they are called by \\ %%wv_info = WV_FN_<name>( [Order, Scale, N ][, /DOUBLE] [, FREQUENCY=variable] [, /SPATIAL] [, WAVELET=variable])%%
    * the 'mother wavelt' is returned to the variable specified with keyword WAVELET
    * scale is a scalar //input// parameter containing the scale of the function (its width if you want)
    * in wavelet is returned the //fourier transform// of the wavelet function, except when keyword %%/SPATIAL%% ist set
    * //gaussian// is a derivative of the gaussian bell where parameter 'order' is the order of the derivative: 1 is the first derivative, 2 is the second derivative (sometimes called the 'mexican hat wavelet')
    * //morlet// is a sine with damped amplitude parameter 'order' is called 'nondimensioanl frequency parameter' in the comments of the function. It is the ratio between amplitude damping constant and frequency. It is a measure for the number of full oscialltions.

In principle you can construct **your own wavelet function** for the continous wavelet transform by writing a function following the name and call conventions:
  * the filename should be %%WV_FN_<name>.PRO%% where <name> is then used as paramter %%family%% when calling %%WV_CWT%%
  * the source code can stand anywhere in the idl path, e.g. in the directory of your current project. 
  * for a continous wavelet you must provide parameters Order, Scale, N, and named paramters WAVELET, SPATIAL, DOUBLE and FREQUENCY.
  * the function must return a structure containing information about the function - allthough it is not evaluated in %%WV_DWT.PRO%%
  * if flag SPATIAL //is set// it must return the wavelet 
  * if flag SPATIAL //is not set// it must return the fourier transform of the wavelet and the related frequencies
But so far i get only weird results trying to implement a continous Haar wavelet ...

==== examples ====

=== Discrete wavelet transform: ===

This is a lot of code - mainly to reconstruct the values of scale and time for the wavelet coefficients.

<code idl>
; generate wavelet function
; info record for wavelet function
info_rec = WV_FN_HAAR( 1, scaling, wavelet, ioff, joff)

; ioff and joff need to be set !
ioff = 0 ; alternativley -N/2+2
joff = 0 ; alternativley -N/2+2

; Compute the discrete wavelet transform
wave_disc = WV_DWT( data, scaling, wavelet,   ioff, joff )

; power spectrum if desired
; wave_disc = sqrt(abs(wave_disc))

; construct scale and time arrays for the wavelet coefficients
di = N/2
; first two coefficients in wv_disc contain the 'residuums' S0 and S1
time_disc = [ 0 , di ] + di/2
scale_disc = [ di , di ]
; the following two coefficients belong to the largest scale
time_disc  = [ time_disc  , time_disc  ]
scale_disc = [ scale_disc , scale_disc ]
; the following coefficents belong to scales 2^i
while di gt 2 do begin
  di /= 2
  for j = 0L, n-1, di do begin
    time_disc = [ time_disc , di/2+j ]
  endfor ; j
  scale_disc = [ scale_disc , replicate(di,n/di) ]
endwhile
; convert to units of time by multiplying with timestep dt
time_disc = time_disc * dt + time[0]
scale_disc *= dt

; make a contour plot 
CONTOUR, $
  wave_disc[2:n-1] , $ ; drop the first two coefficients as they contain the residuums
  time_disc[2:n-1], scale_disc[2:n-1], $ ; reconstructed time and scale from above
   TITLE = 'discrete wavelet transform', $
   XTITLE = 'time', $
   YTITLE = 'period', $  
   /IRREGULAR, $ ; idl will triangulate the data before drawing in a regular grid
   /YLOG , $ ; log scaling on scale axis
   YRANGE = [ max(scale_disc) , min(scale_disc) ], $ ; invert direction of axis
   NLEVELS = 100, $ ; number of isoline levels
   /FILL, $  ; fill isolines with color
   XSTYLE = 1, $
   YSTYLE = 1

</code>

=== Continous wavelet transform ===

<code idl>
; continous wavelet transform
; wave_cont = WV_CWT(data, 'morlet',   8, /PAD, SCALE=scale_cont )

; power spectrum if desired
; wave_cont = sqrt(abs(wave_cont))

; convert to units of time by mutliplying with timestep dt
scale_cont *= dt

; Contour plot
CONTOUR, wave_cont, time, scale_cont, $
   TITLE = 'continous wavelet transform', $
   XTITLE = 'time', $
   YTITLE = 'period', $  
   /YLOG , $ ; log scaling on scale axis
   YRANGE = [ max(scale_cont) , min(scale_cont) ], $ ; invert direction of axis
   NLEVELS = 100, $ ; number of isoline levels
   /FILL, $  ; fill isolines with color
   XSTYLE = 1, $
   YSTYLE = 1

</code>


===== Errors / debugging =====

IDL can handle most of arithmetic errors without stopping. 

You just get a message at the end of the program like:

  % Program caused arithmetic error: Floating divide by 0

or:

  % Program caused arithmetic error: Floating illegal operand

The latter might be something like assigning values outside the int-range to an array which has been defined as   a=intarr(N).

To debug errors you can set the system variable 
  !EXCEPT = 2
this forces IDL to print an error message including the line number after each the line where an error occurs.\\
This should only used for debugging purposes - it slows down IDL by about 5%.

The default is %%!EXCEPT=1%% meaning that error messages are printed at the end of the program. With %%!EXCEPT=0%% no error messages are printed at all.


===== idl code =====

native idl routines and examples can be found in the lib directory of the idl distribution. IT can be found by investigating the [[http://www.physics.nyu.edu/grierlab/idl_html_help/sysvars6.html#wp997094|!PATH]] variable. Currently it is %%/opt/itt/idl8.2/idl82/lib/%%

===== operating system =====

For plotting on the screen and for building pathes your code may need to know which operating system your code is running on.
You can find this out by reading the !VERSION system variable:

on one of our linux machine it looks like:

  IDL> help, !version , /stru
  ** Structure !VERSION, 8 tags, length=104, data length=100:
     ARCH            STRING    'x86_64'
     OS              STRING    'linux'
     OS_FAMILY       STRING    'unix'
     OS_NAME         STRING    'linux'
     RELEASE         STRING    '8.2.1'
     BUILD_DATE      STRING    'Aug 20 2012'
     MEMORY_BITS     INT             64
     FILE_OFFSET_BITS
                     INT             64

On a windows system it may look like
  IDL> Help, !Version, /Structure
   ** Structure !VERSION, 8 tags, length=104, data length=100:
      ARCH            STRING    'x86_64'       OS              STRING    'Win32'
      OS_FAMILY       STRING    'Windows'
      OS_NAME         STRING    'Microsoft Windows'
      RELEASE         STRING    '8.0'
      BUILD_DATE      STRING    'Jun 17 2010'
      MEMORY_BITS     INT             64       FILE_OFFSET_BITS
                      INT             64 

For the correct device To plot on the screen (i.e. open a window to plot in ) you may use:

  if !VERSION.OS_FAMILY eq 'Windows' then set_plot, 'Win' $
                                     else set_plot, 'X'

And similar for pathes:

  if (!version.os_family eq 'linux')  then $
    srtm_dbase = '/home/jschween/exp/ideen/idl/srtm/srtm/' $
  else $ 
    srtm_dbase = 'C:\Users\jschween\Downloads\srtm\'


  
===== IDL version =====

Since November 2012 we have idl 8.2.1

Simply type idl and idl 8.2 will be executed. if no license for idl 8.2 is available 6.4 will be used.
if you want to you use 6.4 explicitly, you need to type idl64. idlhelp and idlde will use the 8.2 versions. for 6.4 use idlhelp64 or idlde64, respectively. 

To directly start version 6.4. type %%/opt/itt/idl6.4/idl64/bin/idl %%

To directly start version 8.2.1. type %%/opt/itt/idl6.4/idl64/bin/idl %%


If your code needs to know under which idl version it is running you may use the !VERSION system variable:

  if !VERSION.RELEASE ge '8.2.1' then begin
    ...
  endif else begin
    ...
  endelse


===== idl license manager =====

To find out how many users are currently in IDL excute:

  /opt/itt/idl8.2/idl82/bin/bin.linux.x86_64/lmstat -c /opt/itt/license/license.dat -f idl

in the listing should appear beside some other ten lines  something like:

  gop: license server UP (MASTER) v11.6
  ...
  Users of idl:  (Total of 60 licenses issued;  Total of 18 licenses in use)  
  ...
   jschween qeb qeb/localhost:11 (v8.2) (gop/1700 522), start Wed 11/11 10:48, 6 licenses
  ...

The line starting with "Users of idl: ..." says that there are 18 license in use. As one allways needs six licences when running idl this means that 18/6=3 (i.e. three) instances of idl are running. The line below indicates that user jschween on machine gop is running one instance of idl (using 6 licenses) v8.2 since Nov. 11.
  
If this is missing and/or you get a long strange output the license manager is probably not running.

To restart the license manager:
  * log in to gop
  * check whether it is running:

  ps -elf | grep idl

  * there should appear the licences manager daemon (lmgrd)
  * if not:

  sudo /etc/init.d/sys5_idl_lmgrd start

password is your password ...

===== dead IDL processes =====

If licencses are obviously blocked by dead idl processes (running for several days, see section about license manager above). You may want to kill them (be careful not to destroy hard work of your colleagues !):
  * Login on the respective machine
  * see what idl processes are running by executing:

  ps -elf | grep idl

  * output may look like:

  ...
  0 S hatpro   18936     1  0  80   0 -  3513 wait    2019 ?        00:00:00 /bin/bash /usr/local/bin/idl /home/dturner/aerioe/tmp/.run_aerioe.pro
  0 R hatpro   18950 18936 98  80   0 - 85711 -       2019 ?        151-17:57:02 /opt/itt/idl8.2/idl82//bin/bin.linux.x86_64/idl /home/dturner/aerioe/tmp/.run_aerioe.pro
  ...

  * indicating that user hatpro has an idl job "%%run_aerioe.pro%%" which was started in 2019 and is running for// 151 days// and ~18hours
  * its process ID is 18950 and it has been started by process 18936
  * eventually check with the (probable) owner whether this is intended 
  * to terminate this process do %%kill 18950%% (this is the process ID in this case !)