text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PersistentVolumeSpec is the specification of a persistent volume.

    IoK8sApiCoreV1PersistentVolumeSpec(;
        accessModes=nothing,
        awsElasticBlockStore=nothing,
        azureDisk=nothing,
        azureFile=nothing,
        capacity=nothing,
        cephfs=nothing,
        cinder=nothing,
        claimRef=nothing,
        csi=nothing,
        fc=nothing,
        flexVolume=nothing,
        flocker=nothing,
        gcePersistentDisk=nothing,
        glusterfs=nothing,
        hostPath=nothing,
        iscsi=nothing,
        __local__=nothing,
        mountOptions=nothing,
        nfs=nothing,
        nodeAffinity=nothing,
        persistentVolumeReclaimPolicy=nothing,
        photonPersistentDisk=nothing,
        portworxVolume=nothing,
        quobyte=nothing,
        rbd=nothing,
        scaleIO=nothing,
        storageClassName=nothing,
        storageos=nothing,
        volumeMode=nothing,
        vsphereVolume=nothing,
    )

    - accessModes::Vector{String} : AccessModes contains all ways the volume can be mounted. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#access-modes
    - awsElasticBlockStore::IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource : AWSElasticBlockStore represents an AWS Disk resource that is attached to a kubelet's host machine and then exposed to the pod. More info: https://kubernetes.io/docs/concepts/storage/volumes#awselasticblockstore
    - azureDisk::IoK8sApiCoreV1AzureDiskVolumeSource : AzureDisk represents an Azure Data Disk mount on the host and bind mount to the pod.
    - azureFile::IoK8sApiCoreV1AzureFilePersistentVolumeSource : AzureFile represents an Azure File Service mount on the host and bind mount to the pod.
    - capacity::Dict{String, IoK8sApimachineryPkgApiResourceQuantity} : A description of the persistent volume's resources and capacity. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#capacity
    - cephfs::IoK8sApiCoreV1CephFSPersistentVolumeSource : CephFS represents a Ceph FS mount on the host that shares a pod's lifetime
    - cinder::IoK8sApiCoreV1CinderPersistentVolumeSource : Cinder represents a cinder volume attached and mounted on kubelets host machine. More info: https://examples.k8s.io/mysql-cinder-pd/README.md
    - claimRef::IoK8sApiCoreV1ObjectReference : ClaimRef is part of a bi-directional binding between PersistentVolume and PersistentVolumeClaim. Expected to be non-nil when bound. claim.VolumeName is the authoritative bind between PV and PVC. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#binding
    - csi::IoK8sApiCoreV1CSIPersistentVolumeSource : CSI represents storage that is handled by an external CSI driver (Beta feature).
    - fc::IoK8sApiCoreV1FCVolumeSource : FC represents a Fibre Channel resource that is attached to a kubelet's host machine and then exposed to the pod.
    - flexVolume::IoK8sApiCoreV1FlexPersistentVolumeSource : FlexVolume represents a generic volume resource that is provisioned/attached using an exec based plugin.
    - flocker::IoK8sApiCoreV1FlockerVolumeSource : Flocker represents a Flocker volume attached to a kubelet's host machine and exposed to the pod for its usage. This depends on the Flocker control service being running
    - gcePersistentDisk::IoK8sApiCoreV1GCEPersistentDiskVolumeSource : GCEPersistentDisk represents a GCE Disk resource that is attached to a kubelet's host machine and then exposed to the pod. Provisioned by an admin. More info: https://kubernetes.io/docs/concepts/storage/volumes#gcepersistentdisk
    - glusterfs::IoK8sApiCoreV1GlusterfsPersistentVolumeSource : Glusterfs represents a Glusterfs volume that is attached to a host and exposed to the pod. Provisioned by an admin. More info: https://examples.k8s.io/volumes/glusterfs/README.md
    - hostPath::IoK8sApiCoreV1HostPathVolumeSource : HostPath represents a directory on the host. Provisioned by a developer or tester. This is useful for single-node development and testing only! On-host storage is not supported in any way and WILL NOT WORK in a multi-node cluster. More info: https://kubernetes.io/docs/concepts/storage/volumes#hostpath
    - iscsi::IoK8sApiCoreV1ISCSIPersistentVolumeSource : ISCSI represents an ISCSI Disk resource that is attached to a kubelet's host machine and then exposed to the pod. Provisioned by an admin.
    - __local__::IoK8sApiCoreV1LocalVolumeSource : Local represents directly-attached storage with node affinity
    - mountOptions::Vector{String} : A list of mount options, e.g. [\"ro\", \"soft\"]. Not validated - mount will simply fail if one is invalid. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes/#mount-options
    - nfs::IoK8sApiCoreV1NFSVolumeSource : NFS represents an NFS mount on the host. Provisioned by an admin. More info: https://kubernetes.io/docs/concepts/storage/volumes#nfs
    - nodeAffinity::IoK8sApiCoreV1VolumeNodeAffinity : NodeAffinity defines constraints that limit what nodes this volume can be accessed from. This field influences the scheduling of pods that use this volume.
    - persistentVolumeReclaimPolicy::String : What happens to a persistent volume when released from its claim. Valid options are Retain (default for manually created PersistentVolumes), Delete (default for dynamically provisioned PersistentVolumes), and Recycle (deprecated). Recycle must be supported by the volume plugin underlying this PersistentVolume. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#reclaiming
    - photonPersistentDisk::IoK8sApiCoreV1PhotonPersistentDiskVolumeSource : PhotonPersistentDisk represents a PhotonController persistent disk attached and mounted on kubelets host machine
    - portworxVolume::IoK8sApiCoreV1PortworxVolumeSource : PortworxVolume represents a portworx volume attached and mounted on kubelets host machine
    - quobyte::IoK8sApiCoreV1QuobyteVolumeSource : Quobyte represents a Quobyte mount on the host that shares a pod's lifetime
    - rbd::IoK8sApiCoreV1RBDPersistentVolumeSource : RBD represents a Rados Block Device mount on the host that shares a pod's lifetime. More info: https://examples.k8s.io/volumes/rbd/README.md
    - scaleIO::IoK8sApiCoreV1ScaleIOPersistentVolumeSource : ScaleIO represents a ScaleIO persistent volume attached and mounted on Kubernetes nodes.
    - storageClassName::String : Name of StorageClass to which this persistent volume belongs. Empty value means that this volume does not belong to any StorageClass.
    - storageos::IoK8sApiCoreV1StorageOSPersistentVolumeSource : StorageOS represents a StorageOS volume that is attached to the kubelet's host machine and mounted into the pod More info: https://examples.k8s.io/volumes/storageos/README.md
    - volumeMode::String : volumeMode defines if a volume is intended to be used with a formatted filesystem or to remain in raw block state. Value of Filesystem is implied when not included in spec. This is a beta feature.
    - vsphereVolume::IoK8sApiCoreV1VsphereVirtualDiskVolumeSource : VsphereVolume represents a vSphere volume attached and mounted on kubelets host machine
"""
mutable struct IoK8sApiCoreV1PersistentVolumeSpec <: SwaggerModel
    accessModes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: accessModes
    awsElasticBlockStore::Any # spec type: Union{ Nothing, IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource } # spec name: awsElasticBlockStore
    azureDisk::Any # spec type: Union{ Nothing, IoK8sApiCoreV1AzureDiskVolumeSource } # spec name: azureDisk
    azureFile::Any # spec type: Union{ Nothing, IoK8sApiCoreV1AzureFilePersistentVolumeSource } # spec name: azureFile
    capacity::Any # spec type: Union{ Nothing, Dict{String, IoK8sApimachineryPkgApiResourceQuantity} } # spec name: capacity
    cephfs::Any # spec type: Union{ Nothing, IoK8sApiCoreV1CephFSPersistentVolumeSource } # spec name: cephfs
    cinder::Any # spec type: Union{ Nothing, IoK8sApiCoreV1CinderPersistentVolumeSource } # spec name: cinder
    claimRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectReference } # spec name: claimRef
    csi::Any # spec type: Union{ Nothing, IoK8sApiCoreV1CSIPersistentVolumeSource } # spec name: csi
    fc::Any # spec type: Union{ Nothing, IoK8sApiCoreV1FCVolumeSource } # spec name: fc
    flexVolume::Any # spec type: Union{ Nothing, IoK8sApiCoreV1FlexPersistentVolumeSource } # spec name: flexVolume
    flocker::Any # spec type: Union{ Nothing, IoK8sApiCoreV1FlockerVolumeSource } # spec name: flocker
    gcePersistentDisk::Any # spec type: Union{ Nothing, IoK8sApiCoreV1GCEPersistentDiskVolumeSource } # spec name: gcePersistentDisk
    glusterfs::Any # spec type: Union{ Nothing, IoK8sApiCoreV1GlusterfsPersistentVolumeSource } # spec name: glusterfs
    hostPath::Any # spec type: Union{ Nothing, IoK8sApiCoreV1HostPathVolumeSource } # spec name: hostPath
    iscsi::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ISCSIPersistentVolumeSource } # spec name: iscsi
    __local__::Any # spec type: Union{ Nothing, IoK8sApiCoreV1LocalVolumeSource } # spec name: local
    mountOptions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: mountOptions
    nfs::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NFSVolumeSource } # spec name: nfs
    nodeAffinity::Any # spec type: Union{ Nothing, IoK8sApiCoreV1VolumeNodeAffinity } # spec name: nodeAffinity
    persistentVolumeReclaimPolicy::Any # spec type: Union{ Nothing, String } # spec name: persistentVolumeReclaimPolicy
    photonPersistentDisk::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PhotonPersistentDiskVolumeSource } # spec name: photonPersistentDisk
    portworxVolume::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PortworxVolumeSource } # spec name: portworxVolume
    quobyte::Any # spec type: Union{ Nothing, IoK8sApiCoreV1QuobyteVolumeSource } # spec name: quobyte
    rbd::Any # spec type: Union{ Nothing, IoK8sApiCoreV1RBDPersistentVolumeSource } # spec name: rbd
    scaleIO::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ScaleIOPersistentVolumeSource } # spec name: scaleIO
    storageClassName::Any # spec type: Union{ Nothing, String } # spec name: storageClassName
    storageos::Any # spec type: Union{ Nothing, IoK8sApiCoreV1StorageOSPersistentVolumeSource } # spec name: storageos
    volumeMode::Any # spec type: Union{ Nothing, String } # spec name: volumeMode
    vsphereVolume::Any # spec type: Union{ Nothing, IoK8sApiCoreV1VsphereVirtualDiskVolumeSource } # spec name: vsphereVolume

    function IoK8sApiCoreV1PersistentVolumeSpec(;accessModes=nothing, awsElasticBlockStore=nothing, azureDisk=nothing, azureFile=nothing, capacity=nothing, cephfs=nothing, cinder=nothing, claimRef=nothing, csi=nothing, fc=nothing, flexVolume=nothing, flocker=nothing, gcePersistentDisk=nothing, glusterfs=nothing, hostPath=nothing, iscsi=nothing, __local__=nothing, mountOptions=nothing, nfs=nothing, nodeAffinity=nothing, persistentVolumeReclaimPolicy=nothing, photonPersistentDisk=nothing, portworxVolume=nothing, quobyte=nothing, rbd=nothing, scaleIO=nothing, storageClassName=nothing, storageos=nothing, volumeMode=nothing, vsphereVolume=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("accessModes"), accessModes)
        setfield!(o, Symbol("accessModes"), accessModes)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("awsElasticBlockStore"), awsElasticBlockStore)
        setfield!(o, Symbol("awsElasticBlockStore"), awsElasticBlockStore)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("azureDisk"), azureDisk)
        setfield!(o, Symbol("azureDisk"), azureDisk)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("azureFile"), azureFile)
        setfield!(o, Symbol("azureFile"), azureFile)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("capacity"), capacity)
        setfield!(o, Symbol("capacity"), capacity)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("cephfs"), cephfs)
        setfield!(o, Symbol("cephfs"), cephfs)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("cinder"), cinder)
        setfield!(o, Symbol("cinder"), cinder)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("claimRef"), claimRef)
        setfield!(o, Symbol("claimRef"), claimRef)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("csi"), csi)
        setfield!(o, Symbol("csi"), csi)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("fc"), fc)
        setfield!(o, Symbol("fc"), fc)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("flexVolume"), flexVolume)
        setfield!(o, Symbol("flexVolume"), flexVolume)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("flocker"), flocker)
        setfield!(o, Symbol("flocker"), flocker)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("gcePersistentDisk"), gcePersistentDisk)
        setfield!(o, Symbol("gcePersistentDisk"), gcePersistentDisk)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("glusterfs"), glusterfs)
        setfield!(o, Symbol("glusterfs"), glusterfs)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("hostPath"), hostPath)
        setfield!(o, Symbol("hostPath"), hostPath)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("iscsi"), iscsi)
        setfield!(o, Symbol("iscsi"), iscsi)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("local"), __local__)
        setfield!(o, Symbol("__local__"), __local__)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("mountOptions"), mountOptions)
        setfield!(o, Symbol("mountOptions"), mountOptions)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("nfs"), nfs)
        setfield!(o, Symbol("nfs"), nfs)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("nodeAffinity"), nodeAffinity)
        setfield!(o, Symbol("nodeAffinity"), nodeAffinity)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("persistentVolumeReclaimPolicy"), persistentVolumeReclaimPolicy)
        setfield!(o, Symbol("persistentVolumeReclaimPolicy"), persistentVolumeReclaimPolicy)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("photonPersistentDisk"), photonPersistentDisk)
        setfield!(o, Symbol("photonPersistentDisk"), photonPersistentDisk)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("portworxVolume"), portworxVolume)
        setfield!(o, Symbol("portworxVolume"), portworxVolume)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("quobyte"), quobyte)
        setfield!(o, Symbol("quobyte"), quobyte)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("rbd"), rbd)
        setfield!(o, Symbol("rbd"), rbd)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("scaleIO"), scaleIO)
        setfield!(o, Symbol("scaleIO"), scaleIO)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("storageClassName"), storageClassName)
        setfield!(o, Symbol("storageClassName"), storageClassName)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("storageos"), storageos)
        setfield!(o, Symbol("storageos"), storageos)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("volumeMode"), volumeMode)
        setfield!(o, Symbol("volumeMode"), volumeMode)
        validate_property(IoK8sApiCoreV1PersistentVolumeSpec, Symbol("vsphereVolume"), vsphereVolume)
        setfield!(o, Symbol("vsphereVolume"), vsphereVolume)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeSpec

const _property_map_IoK8sApiCoreV1PersistentVolumeSpec = Dict{Symbol,Symbol}(Symbol("accessModes")=>Symbol("accessModes"), Symbol("awsElasticBlockStore")=>Symbol("awsElasticBlockStore"), Symbol("azureDisk")=>Symbol("azureDisk"), Symbol("azureFile")=>Symbol("azureFile"), Symbol("capacity")=>Symbol("capacity"), Symbol("cephfs")=>Symbol("cephfs"), Symbol("cinder")=>Symbol("cinder"), Symbol("claimRef")=>Symbol("claimRef"), Symbol("csi")=>Symbol("csi"), Symbol("fc")=>Symbol("fc"), Symbol("flexVolume")=>Symbol("flexVolume"), Symbol("flocker")=>Symbol("flocker"), Symbol("gcePersistentDisk")=>Symbol("gcePersistentDisk"), Symbol("glusterfs")=>Symbol("glusterfs"), Symbol("hostPath")=>Symbol("hostPath"), Symbol("iscsi")=>Symbol("iscsi"), Symbol("local")=>Symbol("__local__"), Symbol("mountOptions")=>Symbol("mountOptions"), Symbol("nfs")=>Symbol("nfs"), Symbol("nodeAffinity")=>Symbol("nodeAffinity"), Symbol("persistentVolumeReclaimPolicy")=>Symbol("persistentVolumeReclaimPolicy"), Symbol("photonPersistentDisk")=>Symbol("photonPersistentDisk"), Symbol("portworxVolume")=>Symbol("portworxVolume"), Symbol("quobyte")=>Symbol("quobyte"), Symbol("rbd")=>Symbol("rbd"), Symbol("scaleIO")=>Symbol("scaleIO"), Symbol("storageClassName")=>Symbol("storageClassName"), Symbol("storageos")=>Symbol("storageos"), Symbol("volumeMode")=>Symbol("volumeMode"), Symbol("vsphereVolume")=>Symbol("vsphereVolume"))
const _property_types_IoK8sApiCoreV1PersistentVolumeSpec = Dict{Symbol,String}(Symbol("accessModes")=>"Vector{String}", Symbol("awsElasticBlockStore")=>"IoK8sApiCoreV1AWSElasticBlockStoreVolumeSource", Symbol("azureDisk")=>"IoK8sApiCoreV1AzureDiskVolumeSource", Symbol("azureFile")=>"IoK8sApiCoreV1AzureFilePersistentVolumeSource", Symbol("capacity")=>"Dict{String, IoK8sApimachineryPkgApiResourceQuantity}", Symbol("cephfs")=>"IoK8sApiCoreV1CephFSPersistentVolumeSource", Symbol("cinder")=>"IoK8sApiCoreV1CinderPersistentVolumeSource", Symbol("claimRef")=>"IoK8sApiCoreV1ObjectReference", Symbol("csi")=>"IoK8sApiCoreV1CSIPersistentVolumeSource", Symbol("fc")=>"IoK8sApiCoreV1FCVolumeSource", Symbol("flexVolume")=>"IoK8sApiCoreV1FlexPersistentVolumeSource", Symbol("flocker")=>"IoK8sApiCoreV1FlockerVolumeSource", Symbol("gcePersistentDisk")=>"IoK8sApiCoreV1GCEPersistentDiskVolumeSource", Symbol("glusterfs")=>"IoK8sApiCoreV1GlusterfsPersistentVolumeSource", Symbol("hostPath")=>"IoK8sApiCoreV1HostPathVolumeSource", Symbol("iscsi")=>"IoK8sApiCoreV1ISCSIPersistentVolumeSource", Symbol("local")=>"IoK8sApiCoreV1LocalVolumeSource", Symbol("mountOptions")=>"Vector{String}", Symbol("nfs")=>"IoK8sApiCoreV1NFSVolumeSource", Symbol("nodeAffinity")=>"IoK8sApiCoreV1VolumeNodeAffinity", Symbol("persistentVolumeReclaimPolicy")=>"String", Symbol("photonPersistentDisk")=>"IoK8sApiCoreV1PhotonPersistentDiskVolumeSource", Symbol("portworxVolume")=>"IoK8sApiCoreV1PortworxVolumeSource", Symbol("quobyte")=>"IoK8sApiCoreV1QuobyteVolumeSource", Symbol("rbd")=>"IoK8sApiCoreV1RBDPersistentVolumeSource", Symbol("scaleIO")=>"IoK8sApiCoreV1ScaleIOPersistentVolumeSource", Symbol("storageClassName")=>"String", Symbol("storageos")=>"IoK8sApiCoreV1StorageOSPersistentVolumeSource", Symbol("volumeMode")=>"String", Symbol("vsphereVolume")=>"IoK8sApiCoreV1VsphereVirtualDiskVolumeSource")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeSpec }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeSpec[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeSpec)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolumeSpec }, name::Symbol, val)
end


text:
module aurora
@inline function run(nt_out, t_trans, D, V, par_loss_rates, src_rad_prof, S_rates,
             R_rates, rr, pro, qpr, r_saw, dlen, time, saw, flx, it_out, dsaw, rcl,
             divbls, taudiv, taupump, tauwret, rvol_lcfs, dbound, dlim, prox,
             rn, alg_opt, evolveneut)
    """
    !
    ! Run forward model of radial impurity transport, returning the density of
    ! each charge state over time and space.
    !
    ! ------ All inputs are in CGS units -----
    !
    ! Get list of required input parameters in Python using
    ! print(aurora.run.__doc__)
    !
    ! Args:
    !     nion         integer
    !                    Number of ionization stages. Not needed for Python calls.
    !     ir           integer
    !                    Number of radial grid points. Not needed for Python calls.
    !     nt           integer
    !                    Number of time steps for the solution. Not needed for Python calls.
    !     nt_out       integer
    !                    Number of times at which the impurity densities shall be saved.
    !     nt_trans     integer
    !                    Number of times at which D,V profiles are given. Not needed for Python calls.
    !     t_trans      real*8 (nt_trans)
    !                    Times at which transport coefficients change [s].
    !     D            real*8 (ir,nt_trans,nion)
    !                    Diffusion coefficient on time and radial grids [cm^2/s]
    !                    This must be given for each charge state and time.
    !     V            real*8 (ir,nt_trans,nion)
    !                    Drift velocity on time and radial grids [cm/s]
    !                    This must be given for each charge state and time.
    !     par_loss_rates  real*8 (ir,nt)
    !                    Frequency for parallel loss on radial and time grids [1/s]
    !     src_rad_prof real*8 (ir,nt)
    !                    Radial profile of neutrals over time: n0(ir,t) = flx*src_rad_prof(ir,t).
    !     S_rates      real*8 (ir,nion,nt)
    !                    Ionisation rates (nz=nion must be filled with zeros).
    !     R_rates      real*8 (ir,nion,nt)
    !                    Recombination rates (nz=nion must be filled with zeros)
    !     rr           real*8 (ir)
    !                    Radial grid, defined using normalized flux surface volumes
    !     pro          real*8 (ir)
    !                    Normalized first derivative of the radial grid, defined by
    !                    pro = (drho/dr)/(2 d_rho) = rho'/(2 d_rho)
    !     qpr          real*8 (ir)
    !                    Normalized second derivative of the radial grid, defined as
    !                    qpr = (d^2 rho/dr^2)/(2 d_rho) = rho''/(2 d_rho)
    !     r_saw        real*8
    !                    Sawteeth inversion radius [cm]
    !     dlen         real*8
    !                    Decay length at last radial grid point
    !     time         real*8 (nt)
    !                    Time grid for transport solver
    !     saw          integer (nt)
    !                    Switch to induce a sawtooth crashes
    !                    If saw(it) eq 1 there is a crash at time(it)
    !     flx          real*8 (nt)
    !                    Impurity influx per unit length and time [1/cm/s]
    !                    (only 1 to nt-1 are used)
    !     it_out       integer (nt)
    !                    Store the impurity distributions if it_out(it).eq.1
    !     dsaw         real*8
    !                    Width of sawtooth crash region.
    !     rcl          real*8
    !                    Wall recycling coefficient. Normally, this would be in the range [0,1].
    !                    However, if set to a value <0, then this is interpreted as a flag, indicating
    !                    that particles in the divertor should NEVER return to the main plasma.
    !     divbls       real*8
    !                    Fraction of the total flux (given by flx) that is directed to the divertor
    !                    reservoir. This value should be between 0 and 1.
    !     taudiv       real*8
    !                    Time scale for transport out of the divertor reservoir [s]
    !     taupump      real*8
    !                    Time scale for impurity elimination through out-pumping [s]
    !     tauwret      real*8
    !                    Time scale of temporary retention at the wall [s]
    !     rvol_lcfs    real*8
    !                    Radius (in rvol units, cm) at which the LCFS is located
    !     dbound       real*8
    !                    Width of the SOL, given by r_bound - r_lcfs (in rvol coordinates, cm)
    !                    This value sets the width of the radial grid.
    !     dlim         real*8
    !                    Position of the limiter wrt to the LCFS, i.e. r_lim - r_lcfs (cm, in r_vol units).
    !                    Inside of this limiter location, the parallel connection length to the divertor applies,
    !                    while outside of it the relevant connection length is the one to the limiter.
    !                    These different connection lengths must be taken into consideration when
    !                    preparing the parallel loss rate variable.
    !     prox         real*8
    !                    Grid parameter for loss rate at the last radial point, returned by
    !                    `get_radial_grid' subroutine.
    !     rn_t0        real*8 (ir,nion), optional
    !                    Impurity densities at the start time [1/cm^3]. If not provided, all elements are
    !                    set to 0.
    !     alg_opt      integer, optional
    !                    Integer to indicate algorithm to be used.
    !                    If set to 0, use the finite-differences algorithm used in the 2018 version of STRAHL.
    !                    If set to 1, use the Linder finite-volume algorithm (see Linder et al. NF 2020)
    !     evolneut     logical, optional
    !                    Boolean to activate evolution of neutrals (like any ionization stage)
    !
    ! Returns:
    !
    !     rn_out       real*8 (ir,nion,nt_out)
    !                    Impurity densities (temporarily) in the magnetically-confined plasma at the
    !                    requested times [1/cm^3].
    !     N_ret        real*8 (nt_out)
    !                    Impurity densities (permanently) retained at the wall over time [1/cm^3].
    !     N_wall       real*8 (nt_out)
    !                    Impurity densities (temporarily) at the wall over time [1/cm^3].
    !     N_div        real*8 (nt_out)
    !                    Impurity densities (temporarily) in the divertor reservoir over time [1/cm^3].
    !     N_pump       real*8 (nt_out)
    !                    Impurity densities (permanently) in the pump [1/cm^3].
    !     N_tsu        real*8 (nt_out)
    !                    Edge loss [1/cm^3].
    !     N_dsu        real*8 (nt_out)
    !                    Parallel loss [1/cm^3].
    !     N_dsul       real*8 (nt_out)
    !                    Parallel loss to limiter [1/cm^3].
    !     rcld_rate    real*8 (nt_out)
    !                    Recycling from divertor [1/cm^3/s].
    !     rclw_rate    real*8 (nt_out)
    !                    Recycling from wall [1/cm^3/s].
    ! ---------------------------------------------------------------------------

    """
    # Sizing info
    ir, nt_trans, nion = size(D)
    _, nt = size(par_loss_rates)

    # Intermediaries
    ra = Array{Float64,2}(undef, ir, nion)
    diff = Array{Float64,2}(undef, ir, nion)
    conv = Array{Float64,2}(undef, ir, nion)
    sext = Vector{Float64}(undef, ir)

    if alg_opt == 1
        # For use in impden_conserv!
        ntt = Vector{Float64}(undef, ir)
        a = Vector{Float64}(undef, ir)
        b = Vector{Float64}(undef, ir)
        c = Vector{Float64}(undef, ir)
        d = Vector{Float64}(undef, ir)
    else
        # For use in impden!
        a = Array{Float64,2}(undef, ir, nion)
        b = Array{Float64,2}(undef, ir, nion)
        c = Array{Float64,2}(undef, ir, nion)
        d1 = Vector{Float64}(undef, ir)
        bet = Vector{Float64}(undef, ir)
        gam = Vector{Float64}(undef, ir)
    end

    # Outputs
    rn_out = Array{Float64,3}(undef, ir, nion, nt_out)
    N_wall = Vector{Float64}(undef, nt_out)
    N_div = Vector{Float64}(undef, nt_out)
    N_pump = Vector{Float64}(undef, nt_out)
    N_ret = Vector{Float64}(undef, nt_out)
    N_tsu = Vector{Float64}(undef, nt_out)
    N_dsu = Vector{Float64}(undef, nt_out)
    N_dsul = Vector{Float64}(undef, nt_out)
    rcld_rate = Vector{Float64}(undef, nt_out)
    rclw_rate = Vector{Float64}(undef, nt_out)

    rcld = 0.
    rclw = 0.
    Nret = 0.
    tve = 0.
    divnew = 0.
    npump = 0.
    tsu = 0.
    dsu = 0.
    dsul = 0.

    it = 1
    kt = 1
    # Output if indicated
    @inbounds if it_out[it] == 1
        rn_out[:,:,kt] = rn

        N_wall[kt] = tve
        N_div[kt] = divnew
        N_pump[kt] = npump
        N_tsu[kt] = tsu
        N_dsu[kt] = dsu
        N_dsul[kt] = dsul

        N_ret[kt] = Nret
        rcld_rate[kt] = 0.
        rclw_rate[kt] = 0.
        kt += 1
    end

    # ============ time loop: ==========
    @inbounds for it in 2:nt
        dt = time[it] - time[it-1]

        ra .= rn # update old array to new (from previous time step)

        for nz in 1:nion
            # updated transport coefficients for each charge state
            linip_arr!(nt_trans, ir, t_trans, view(D, :, :, nz), time[it], view(diff, :, nz))
            linip_arr!(nt_trans, ir, t_trans, view(V, :, :, nz), time[it], view(conv, :, nz))
        end

        divold = divnew

        # pick current ioniz and recomb coeffs. || loss rate and source radial prof
        r_rates = view(R_rates,:,:,it)
        s_rates = view(S_rates,:,:,it)
        par_loss_rate = view(par_loss_rates,:,it)
        src_prof = view(src_rad_prof,:,it)

        # evolve impurity density with current transport coeffs
        if alg_opt == 1
            Nret, rcld, rclw = impden_conserv!(nion, ir, ra, rn, diff, conv,
                               par_loss_rate, src_prof, s_rates, r_rates, rr,
                               flx[it-1], dlen, dt, rcl, tsu, dsul, divold,
                               divbls, taudiv, tauwret, evolveneut, Nret, rcld,
                               rclw, ntt, a, b, c, d)
        else
            Nret, rcld, rclw = impden!(nion, ir, ra, rn, diff, conv, par_loss_rate, src_prof, s_rates, r_rates, rr,
                            pro, qpr, flx[it-1], dlen, dt, rcl, tsu, dsul, divold,
                            divbls, taudiv, tauwret, Nret, rcld, rclw, a, b, c, d1, bet, gam)
        end

        # sawteeth
        if saw[it] == 1
            saw_mix!(nion, ir, rn, r_saw, dsaw, rr, pro)
        end

        # particle losses at wall & divertor + recycling
        tve, divnew, npump, tsu, dsu, dsul = edge_model!(nion, ir, ra, rn, diff, conv, par_loss_rate, dt, rvol_lcfs,
                            dbound, dlim, prox, rr, pro, rcl, taudiv, taupump,
                            divbls, divold, flx[it-1], divnew, tve, npump, tsu,
                            dsu, dsul)

        # array time-step saving/output
        if it_out[it] == 1
            rn_out[:,:,kt] = rn

            N_wall[kt] = tve
            N_div[kt] = divnew
            N_pump[kt] = npump
            N_tsu[kt] = tsu
            N_dsu[kt] = dsu
            N_dsul[kt] = dsul

            N_ret[kt] = Nret
            rcld_rate[kt] = rcld
            rclw_rate[kt] = rclw
            kt += 1
        end
    end

    # End of time loop

    return rn_out, N_wall, N_div, N_pump, N_ret, N_tsu, N_dsu, N_dsul, rcld_rate, rclw_rate
end
@inline function impden_conserv!(nion, ir, ra, rn, diff, conv, par_loss_rate, src_prof, s_rates, r_rates, rr,
                flx, fall_outsol, det, rcl, tsuold, dsulold, divold, divbls, taudiv,
                tauwret, evolveneut, Nret, rcld, rclw, nt, a, b, c, d)
    #     **************************************************
    #     ** Time CENTERED TREATMENT in radial coordinate **
    #     ** LACKNER METHOD IN Z                          **
    #     ** time step is split in two half time steps    **
    #     **************************************************
    flxtot = flx*(1-divbls)
    dt = det/2.0

    # Recycling
    if rcl >= 0    # activated divertor return (R>=0) + recycling mode (if R>0)
        rcld = divold/taudiv
        rclw = rcl*(tsuold+dsulold)

        if tauwret > 0
          # include flux from particles previously retained at the wall
          Nret = Nret * (1-det/tauwret) + rclw*det     # number of particles temporarily retained at wall
          rclw = Nret/tauwret    # component that goes back to be a source
        end

        flx_rcl = rclw + rcld
    else   # no divertor return at all
        rcld = 0.0
        rclw = 0.0
        flx_rcl = 0.0
    end

    # Radial profile of recycling source should be close to the wall
    # Set it to be the same as the input (valve) source -- this may need revision:
    srcl = src_prof

    # select whether neutrals should be evolved
    ns = 2
    if evolveneut
        ns = 1
    end

    # ----- First half time step: direction up --------------------|
    # Externally provided neutrals
    if ns == 1
       impden_constTranspMatrix(view(ra,:,1), view(diff,:,1), view(conv,:,1), ir, dt, fall_outsol, par_loss_rate, rr, a, b, c, d)

       for i=1:ir
          b[i]    = b[i] + dt*s_rates[i,1]
          d[i]    = d[i] + dt*(flxtot*src_prof[i] + flx_rcl*srcl[i])
       end
       TDMA(a, b, c, d, ir, nt)
    else
       nt .= 0
       # radial profile of neutrals (given as input)
       for i=1:ir
           rn[i,1] = flxtot*src_prof[i] + flx_rcl* srcl[i]
       end
    end

    # Ions and recombined neutrals
    for nz=ns:nion
       # Construct transport matrix
       impden_constTranspMatrix(view(ra,:,nz), view(diff,:,nz), view(conv,:,nz), ir, dt, fall_outsol, par_loss_rate, rr, a, b, c, d)

       #Add ionization and recombination
       for i=1:ir
          b[i]    = b[i] + dt*s_rates[i,nz]
          if nz > 1
               d[i]    = d[i] + dt*(rn[i,nz-1]*s_rates[i,nz-1] - ra[i,nz]*r_rates[i,nz-1])
          end
          if nz < nion
               d[i]    = d[i] + dt*ra[i,nz+1]*r_rates[i,nz]
          end
          if nz == 2
               d[i]    = d[i] + dt*nt[i]*s_rates[i,1]
          end
       end

       # Solve tridiagonal system of equations
       TDMA(a, b, c, d, ir, view(rn,:,nz))
    end

    # ----- Second half time step: direction down -----------------|
    # Ions and recombined neutrals
    for nz=nion:-1:ns
       # Construct transport matrix
       impden_constTranspMatrix(view(rn,:,nz), view(diff,:,nz), view(conv,:,nz), ir, dt, fall_outsol, par_loss_rate, rr, a, b, c, d)

       # Add ionization and recombination
       for i=1:ir
          d[i]    = d[i] - dt*rn[i,nz]*s_rates[i,nz]
          if nz > 1
             b[i]    = b[i] + dt*r_rates[i,nz-1]
             d[i]    = d[i] + dt*rn[i,nz-1]*s_rates[i,nz-1]
          end
          if nz < nion
               d[i]    = d[i] + dt*rn[i,nz+1]*r_rates[i,nz]
          end
          if nz == 2
               d[i]    = d[i] + dt*nt[i]*s_rates[i,1]
          end
       end

       # Solve tridiagonal equation system
       TDMA(a, b, c, d, ir, view(rn,:,nz))
    end

    # Externally provided neutrals
    if ns == 1
       impden_constTranspMatrix(nt, view(diff,:,1), view(conv,:,1), ir, dt, fall_outsol, par_loss_rate, rr, a, b, c, d)

       for i=1:ir
          d[i]    = d[i] - dt*nt[i]*s_rates[i,1] + dt*(flxtot*src_prof[i] + flx_rcl* srcl[i])
       end
       TDMA(a, b, c, d, ir, view(rn,:,1))
    end

    return Nret, rcld, rclw
end

@inline function impden_constTranspMatrix(rnt, dimp, vimp, ir, dt, fall_outsol, par_loss_rate, rr, a, b, c, d)

    # Temporal component
    a .= 0.0
    b .= 1.0
    c .= 0.0
    d .= rnt

    # Central point (r = 0)
    # Enforces dn/dr|_r=0 = 0 and v(r=0) = 0

    # Geometric factors
    drp = rr[2] - rr[1]
    rp = 0.5*(rr[1] + rr[2])
    gi = 2.0/(rp^2)

    # Diffusion contribution
    Dp = 0.5*(dimp[1] + dimp[2])
    coefp = 0.5*dt*gi*rp*Dp/drp

    b[1] += coefp
    c[1] -= coefp
    d[1] += coefp*(rnt[2] - rnt[1])

    # adeection contribution
    vp = 0.5*(vimp[1]+vimp[2])
    kap = 0.0
    if abs(vp) > 2.0*Dp/drp
        kap = max(0.0, 1.0-2.0*Dp/(drp*abs(vp)))*sign(vp)
    end

    coefp = 0.25*dt*gi*rp*vp

    b[1] += coefp*(1 + kap)
    c[1] += coefp*(1 - kap)
    d[1] -= coefp*((1 + kap)*rnt[1] + (1 - kap)*rnt[2])

    # Interior points
    for i in 2:ir-1
        #! ----- Geometric factors -------------------------------------|
        drm     = drp
        drp     = rr[i+1]-rr[i]
        rm      = rp
        rp      = .5*(rr[i]+rr[i+1])
        gi      = 2.0/(rp^2-rm^2)

        #! ----- Diffusion contribution --------------------------------|
        Dm      = Dp
        Dp      = .5*(dimp[i]+dimp[i+1])
        coefm   = .5*dt*gi*rm*Dm/drm
        coefp   = .5*dt*gi*rp*Dp/drp

        a[i]    = a[i] - coefm
        b[i]    = b[i] + coefm + coefp
        c[i]    = c[i]         - coefp
        d[i]    = d[i] - coefm*(rnt[i] - rnt[i-1]) + coefp*(rnt[i+1] - rnt[i])

        #! ----- adeection contribution --------------------------------|
        vm      = vp
        vp      = .5*(vimp[i]+vimp[i+1])

        kam     = 0.
        if abs(vm) > 2.0*Dm/drm
           kam = sign(vm)*max(0, 1.0-2.0*Dm/(drm*abs(vm)))
        end
        kap     = 0.
        if abs(vp) > 2.0*Dp/drp
           kap = sign(vp)*max(0, 1.0-2.0*Dp/(drp*abs(vp)))
        end

        coefm   = .25*dt*gi*rm*vm
        coefp   = .25*dt*gi*rp*vp

        a[i]    = a[i] - coefm*(1+kam)
        b[i]    = b[i] - coefm*(1-kam) + coefp*(1+kap)
        c[i]    = c[i]                 + coefp*(1-kap)
        d[i]    = d[i] + coefm*((1+kam)*rnt[i-1] + (1-kam)*rnt[i])- coefp*((1+kap)*rnt[i] + (1-kap)*rnt[i+1])
    end

    # Outer point (r = r)

    # Geometric factors
    drm = drp
    rm = rp
    rp = rr[ir] + 0.5*drm
    gi = 2.0/(rp^2-rm^2)

    # Diffusion contribution
    Dm = Dp
    Dp = dimp[ir]
    coefm = 0.5*dt*gi*rm*Dm/drm
    coefp = 0.5*dt*gi*rp*Dp

    a[ir] -= coefm
    b[ir] += coefm
    d[ir] -= coefm*(rnt[ir] - rnt[ir-1])

    if fall_outsol > 0.0
        b[ir] += coefp/fall_outsol
        d[ir] -= coefp/fall_outsol*rnt[ir]
    end

    # adeection contribution
    vm = vp
    vp = max(0.0, vimp[ir])

    kam = 0.0
    if abs(vm) > 2.0*Dm/drm
        kam = max(0.0, 1.0-2.0*Dm/(drm*abs(vm)))*sign(vm)
    end
    kap = 0.0
    if abs(vp) > 2.0*Dp/drp
        kap = max(0.0, 1.0-2.0*Dp/(drp*abs(vp)))*sign(vp)
    end
    coefm = .25*dt*gi*rm*vm
    coefp = .25*dt*gi*rp*vp

    a[ir] -= coefm*(1+kam)
    b[ir] -= coefm*(1-kam)
    d[ir] += coefm*((1+kam)*rnt[ir-1] + (1-kam)*rnt[ir])

    if fall_outsol > 0
        b[ir] += coefp*(2-(1-kap)*drp/fall_outsol)
        d[ir] -= coefp*(2-(1-kap)*drp/fall_outsol)*rnt[ir]
    end

    # SOL LOSSES

    # Interior points
    for i in 2:ir-1
        coefm = 0.5*dt*par_loss_rate[i]
        b[i] += coefm
        d[i] -= coefm*rnt[i]
    end

    # Outer point (r = r)
    if fall_outsol > 0.0
        coefm   = .5*dt*par_loss_rate[ir]
        b[ir] += coefm
        d[ir] -= coefm*rnt[ir]
   end
end

@inline function TDMA(dLow, dMain, dUpp, rhs, n, sol)
    dMainP = similar(dMain)
    rhsP = similar(rhs)
    # Modify coefficients
    dMainP[1] = dMain[1]
    rhsP[1] = rhs[1]
    for i in 2:n
        dMainP[i] = dMain[i] - dLow[i]*dUpp[i-1]/dMainP[i-1]
        rhsP[i] = rhs[i] - dLow[i]*rhsP[i-1]/dMainP[i-1]
    end

    # Construct solution
    sol[n] = rhsP[n]/dMainP[n]
    for i in n-1:-1:1
        sol[i] = (rhsP[i] - dUpp[i]*sol[i+1])/dMainP[i]
    end
end
@inline function impden!(nion, ir, ra, rn, diff, conv, par_loss_rate, src_prof, s_rates, r_rates, rr, pro, qpr,
                flx, dlen, det, rcl, tsuold, dsulold, divold, divbls, taudiv,
                tauwret, Nret, rcld, rclw, a, b, c, d1, bet, gam)
    """

    """

    # dr near magnetic axis

    der = rr[2] - rr[1]

    #     **************************************************
    #     ** Time CENTERED TREATMENT in radial coordinate **
    #     ** LACKNER METHOD IN Z                          **
    #     ** time step is split in two half time steps    **
    #     **************************************************
    flxtot = flx

    # ------ Recycling ------  ! impden.f  L328-L345
    # Part of the particles that hit the wall are taken to be fully-stuck. Another part is only temporarily retained at the wall.
    # Particles FULLY STUCK to the wall will never leave, i.e. tve can only increase over time (see above).
    # Particles that are only temporarily retained at the wall (given by the recycling fraction) come back
    # (recycle) according to the tauwret time scale.

    if rcl >= 0. # activated divertor return (R>=0) + recycling mode (if R>0)
        rcld = divold/taudiv
        rclw = rcl*(tsuold+dsulold)

        if tauwret > 0.
            # include flux from particles previously retained at the wall
            Nret = Nret * (1 - det/tauwret) + rclw*det # number of particles temporarily retained at wall
            rclw = Nret/tauwret # component that goes back to be a source
        end

        flxtot = flx*(1-divbls) + rclw + rcld
    else # no divertor return at all
        rcld = 0.
        rclw = 0.
    end

    @inbounds for i in 1:ir
        rn[i,1] = flxtot*src_prof[i] # index of 1 stands for first ionization stage
    end

    dt = det/2.

    ########### First half time step direction up #########

    @inbounds for nz in 2:nion
        #   r = 0
        a[1,nz] = 0.
        c[1,nz] = -2.0*dt*diff[1,nz]*pro[1]
        b[1,nz] = 1.0 - c[1,nz] + 2.0*dt*conv[2,nz]/der
        d1[1] = ra[1,nz]*(2.0-b[1,nz])-ra[2,nz]*c[1,nz]

        #   r=r or at r+db respectively

        if dlen > 0.
            temp1 = 4.0*dt*pro[ir]^2*diff[ir,nz]
            temp2 = 0.5*dt*(qpr[ir]*diff[ir,nz]-pro[ir]*conv[ir,nz])
            temp3 = 0.5*dt*(par_loss_rate[ir] + conv[ir,nz]/rr[ir])
            temp4 = 1.0/pro[ir]/dlen

            a[ir,nz] = -temp1
            b[ir,nz] = 1.0+(1.0+0.5*temp4)*temp1 + temp4*temp2 + temp3
            c[ir,nz] = 0.
            d1[ir] = -ra[ir-1,nz]*a[ir,nz] + ra[ir,nz]*(2.0 - b[ir,nz])
            b[ir,nz] += dt*s_rates[ir,nz]
            d1[ir] -= dt*(ra[ir,nz]*r_rates[ir,nz-1] - rn[ir,nz-1]*s_rates[ir,nz-1])

            if nz < nion
                d1[ir] += dt*ra[ir,nz+1]*r_rates[ir,nz]
            end
        else # Edge Conditions for rn[ir]=0
            a[ir,nz] = 0
            b[ir,nz] = 1.
            c[ir,nz] = 0.
            d1[ir] = ra[ir,nz]
        end

        #   normal coefficients

        for i in 2:ir-1
            temp1 = dt*pro[i]^2
            temp2 = 4.0*temp1*diff[i,nz]
            temp3 = dt/2.0*(par_loss_rate[i] + pro[i]*(conv[i+1,nz] - conv[i-1, nz]) + conv[i,nz]/rr[i])
            a[i,nz] = (0.5*dt*qpr[i]*diff[i,nz] + temp1*(0.5*(diff[i+1,nz] -
                     diff[i-1,nz] - conv[i,nz]/pro[i]) - 2.0*diff[i,nz]))
            b[i,nz] = 1.0 + temp2 + temp3
            c[i,nz] = -temp2 - a[i,nz]
            d1[i] = -ra[i-1,nz]*a[i,nz] + ra[i,nz]*(2.0-b[i,nz]) - ra[i+1,nz]*c[i,nz]
        end

        for i in 1:ir-1
            b[i,nz] += dt*s_rates[i,nz]
            d1[i] -= dt*(ra[i,nz]*r_rates[i,nz-1] - rn[i,nz-1]*s_rates[i,nz-1]) # ra --> "alt" (i-1) density; rnz --> "current" [i] density

            if nz < nion
                d1[i] += dt*ra[i,nz+1]*r_rates[i,nz] # al --> recomb; s --> ioniz
            end
        end

        #       solution of tridiagonal equation system
        bet[1] = b[1,nz]
        gam[1] = d1[1]/b[1,nz]
        for i in 2:ir
            bet[i] = b[i,nz] - (a[i,nz]*c[i-1,nz])/bet[i-1]
            gam[i] = (d1[i] - a[i,nz]*gam[i-1])/bet[i]
        end
        rn[ir,nz] = gam[ir]

        for i in ir-1:-1:1
            rn[i,nz] = gam[i] - (c[i,nz]*rn[i+1,nz])/bet[i]
        end
    end

    ########### Second half time step direction down #########

    @inbounds for nz in nion:-1:2
        #   r = 0
        a[1,nz] = 0.
        c[1,nz] = -2.0*dt*diff[1,nz]*pro[1]
        b[1,nz] = 1.0 - c[1,nz] + 2.0*dt*conv[2,nz]/der
        d1[1] = rn[1,nz]*(2.0 - b[1,nz]) - rn[2,nz]*c[1,nz]

        #       r=r or at r+db respectively

        if dlen > 0.
            temp1 = 4.0*dt*pro[ir]^2*diff[ir,nz]
            temp2 = 0.5*dt*(qpr[ir]*diff[ir,nz] - pro[ir]*conv[ir,nz])
            temp3 = 0.5*dt*(par_loss_rate[ir] + conv[ir,nz]/rr[ir])
            temp4 = 1.0/pro[ir]/dlen

            a[ir,nz] = -temp1
            b[ir,nz] = 1.0 + (1.0 + 0.5*temp4)*temp1 + temp4*temp2 + temp3
            c[ir,nz] = 0.
            d1[ir] = -rn[ir-1,nz]*a[ir,nz] + rn[ir,nz]*(2.0 - b[ir,nz])
            b[ir,nz] += dt*r_rates[ir,nz-1]
            d1[ir] -= dt*(rn[ir,nz]*s_rates[ir,nz] - rn[ir,nz-1]*s_rates[ir,nz-1])

            if nz < nion
                d1[ir] += dt*rn[ir,nz+1]*r_rates[ir,nz]
            end
        else # Edge conditions for rn[ir]=0.
            a[ir,nz] = 0
            b[ir,nz] = 1.
            c[ir,nz] = 0.
            d1[ir] = rn[ir,nz]
        end

        #       normal coefficients

        for i in 2:ir-1
            temp1 = dt*pro[i]^2
            temp2 = 4.0*temp1*diff[i,nz]
            temp3 = dt/2.0*(par_loss_rate[i] + pro[i]*(conv[i+1,nz] - conv[i-1,nz]) + conv[i,nz]/rr[i])

            a[i,nz] = (0.5*dt*qpr[i]*diff[i,nz] + temp1*(0.5*(diff[i+1,nz]-
                      diff[i-1,nz] - conv[i,nz]/pro[i]) - 2.0*diff[i,nz]))
            b[i,nz] = 1.0 + temp2 + temp3
            c[i,nz] = -temp2 - a[i,nz]
            d1[i] = -rn[i-1,nz]*a[i,nz] + rn[i,nz]*(2.0 - b[i,nz]) - rn[i+1,nz]*c[i,nz]
        end

        for i in 1:ir-1
            b[i,nz] += dt*r_rates[i,nz-1]
            d1[i] -= dt*(rn[i,nz]*s_rates[i,nz] - rn[i,nz-1]*s_rates[i,nz-1])

            if nz < nion
                d1[i] += dt*rn[i,nz+1]*r_rates[i,nz]*(nz < nion)
            end
        end

        #   solution of tridiagonal equation system

        bet[1] = b[1,nz]
        gam[1] = d1[1]/b[1,nz]
        for i in 2:ir
            bet[i] = b[i,nz] - (a[i,nz]*c[i-1,nz])/bet[i-1]
            gam[i] = (d1[i] - a[i,nz]*gam[i-1])/bet[i]
        end
        rn[ir,nz] = gam[ir]

        for i in ir-1:-1:1
            rn[i,nz] = gam[i] - (c[i,nz]*rn[i+1,nz])/bet[i]
        end
    end
    return Nret, rcld, rclw
end

@inline function saw_mix!(nion, ir, rn, rsaw, dsaw, rr, pro)
    """

    """
    #   index of mixing radius
    @inbounds for imix in 1:ir
        if rr[i] > rsaw

        end
    end

    @inbounds for nz in 2:nion    # loop over ionized stages

        #   area integral in mixing radius of old profile

        sum_old = (0.125*(rn[imix,nz]*rr[imix]/pro[imix] -
                   rn[imix-1,nz]*rr[imix-1]/pro[imix-1])) # only use new density, rn
        for i in 2:imix-1
            sum_old += rn[i,nz]*rr[i]/pro[i]
        end

        #   ERFC sawtooth crash model
        ff = sum_old/rr[imix]^2 # nmean
        for i in 1:ir
            rn[i,nz] = (ff/2.0*erfc((rr[i] - rsaw)/dsaw) + (rn[i,nz]/2.0)*
                        erfc(-(rr[i]-rsaw)/dsaw))
        end

        #      flat profile
        #  ff = sum_old/rr(imix)**2
        #  do i=1,imix-1
        #    rn[i,nz] = ff
        #  end do
        #  rn(imix,nz) = (ra(imix+1,nz)+ff)/2.

        #      area integral in mixing radius of new profile
        sum = (0.125*(rn[imix,nz]*rr[imix]/pro[imix] -
               rn[imix-1,nz]*rr[imix-1]/pro[imix-1]))
        for i in 2:imix-1
            sum += rn[i,nz]*rr[i]/pro[i]
        end

        #   ensure particle conservation

        ff = sum_old/sum
        for i in 1:imix
            rn[i,nz] = rn[i,nz]*ff
        end
    end
end

@inline function edge_model!(nion, ir, ra, rn, diff, conv, par_loss_rate, det, rvol_lcfs, dbound,
                    dlim, prox, rr, pro, rcl, taudiv, taupump, divbls, divold,
                    flx, divnew, tve, npump, tsu, dsu, dsul)
    """

    """

    ids = 0
    idl = 0
    rx = rvol_lcfs + dbound     # wall (final) location

    # --------------------------------
    # from plad.f
    @inbounds for i in 2:ir
        if rr[i] <= rvol_lcfs
            ids = i + 1 # number of radial points inside of LCFS
        end
        if rr[i] <= (rvol_lcfs + dlim)
            idl = i + 1 # number of radial points inside of limiter
        end
    end

    # ---------------------------------------
    # ions lost at periphery (not parallel)
    tsu = 0.
    @inbounds for nz in 2:nion
        tsu -= (prox*(diff[ir-1,nz] + diff[ir,nz])*(rn[ir,nz] + ra[ir,nz] -
                rn[ir-1,nz] - ra[ir-1,nz]) + 0.5*(conv[ir-1,nz] + conv[ir,nz])*
                (rn[ir,nz] + ra[ir,nz] + rn[ir-1,nz] + ra[ir-1,nz]))
    end
    tsu = tsu*0.5*pi*rx

    # parallel losses / second

    dsu = 0.
    dsul = 0.

    @inbounds for nz in 2:nion
        for i in ids:idl-1
            dsu += (ra[i,nz] + rn[i,nz])*par_loss_rate[i]*rr[i]/pro[i]
        end
        for i in idl:ir-1
            dsul += (ra[i,nz] + rn[i,nz])*par_loss_rate[i]*rr[i]/pro[i]
        end
    end
    dsu = dsu*pi/2. # to divertor
    dsul = dsul*pi/2.   # to limiter

    # time integrated losses at wall/limiters
    tve += (dsul + tsu)*(1.0 - max(0, rcl))*det

    # particles in divertor
    # If recycling is on, particles from limiter and wall come back.
    # Particles in divertor can only return (with rate given by N_div/taudiv) if rcl>=0
    if rcl >= 0 # activated divertor return (rcl>=0) + recycling mode (if rcl>0)
        taustar = 1.0/(1.0/taudiv + 1.0/taupump)
        ff = 0.5*det/taustar
        divnew = (divold*(1.0-ff) + (dsu + flx*divbls)*det)/(1.0+ff)    # FS corrected
    else
        divnew = divold + (dsu + flx*divbls)*det
    end

    # particles in pump
    npump += 0.5*(divnew + divold)/taupump*det
    return tve, divnew, npump, tsu, dsu, dsul
end

@inline function linip_arr!(m, n, ts, ys, t, y)
    """
    !
    !     Linear interpolation of time dependent profiles.
    !     For times which are outside the bounds the profile is set to the last/first present time.
    !

    !
    !     INPUTS:
    !     m                       number of ts-vector elements
    !     n                        number of used rho-vector elements
    !     ts	                The time vector. Values MUST be monotonically increasing.
    !     ys(rho, ts)          array of ordinate values.
    !     t	                The time for which the ordinates shall be interpolated
    !
    !     OUTPUT:
    !     y(rho)                 interpolated profiles at the chosen time
    !

    """

    #if only one time point has been given, then return the same/only value
    #or if the requested time is greater than the maximum time provided, then use the last time
     @inbounds if m == 1 || t >= ts[m]
        for i in 1:n
           y[i] = ys[i,m]
        end
        return
     end

     #if the requested time is smaller than the minimum time provided, then use the first time
     @inbounds if t < ts[1]
        for i in 1:n
           y[i] = ys[i,1]
        end
        return
     end

     #find time point in "ts" vector that is marginally larger than the requested "t"
     @inbounds for k in 1:m+1
        if (t < ts[k+1])

        end
     end

     ff = (t - ts[k])/(ts[k+1] - ts[k])
     @inbounds for i in 1:n
        y[i] = ys[i,k]*(1.0 - ff) + ys[i,k+1]*ff
     end
end
end

ext: jl

max_stars_repo_name: ToFuProject/Aurora

max_issues_repo_name: ToFuProject/Aurora

max_forks_repo_name: ToFuProject/Aurora


text:
"""
    husimi_Q(rho_in::Array{T1,2},xgrid::Array{Float64,1},
                                      ygrid::Array{Float64,1}) where {T1 <: Number}

Generates the Husimi Q distribution of a quantum state on a rectangular grid in phase space

## args
* rho_in:   d x d array, the density matrix of the state
* xgrid:    vector of floats, x range to calculate the Q function at
* ygrid:    vector of floats, y range to calculate the Q function at

## returns
* Q:        matrix of the calculated Q function
"""
function husimi_Q(rho_in::Array{T1,2},xgrid::Array{Float64,1},
                                      ygrid::Array{Float64,1}) where {T1 <: Number}

    dim = Int(maximum([3*maximum(abs.(xgrid))^2,3*maximum(abs.(ygrid))^2,size(rho_in,1)]))

    # Disp = zeros(size(rho_in))
    Q = zeros(size(xgrid,1),size(ygrid,1))

    lower, raise, N = create_HO(dim)

    rho_temp = zeros(dim,dim)
    rho_temp = convert(Array{ComplexF64,2},rho_temp)
    rho_temp[1:size(rho_in,1),1:size(rho_in,1)] = rho_in

    state = zeros(dim,1)
    state[1] = 1

    for ii = 1:1:size(xgrid,1)
        for jj = 1:1:size(ygrid,1)
            @inbounds temp_a = xgrid[ii]+1im*ygrid[jj]
            Disp = exp(temp_a*raise - conj(temp_a)*lower)
            temp_v = Disp*state
            @inbounds Q[ii,jj] = real(temp_v'*rho_temp*temp_v./pi)[1]
            #Q[ii:ii,jj:jj] = real(temp_v'*rho_temp*temp_v./pi) also works
        end
    end

    return Q
end

"""
    husimi_Q_pw(rho_in::Array{T1,2},
                xgrid::Array{Float64,1},
                ygrid::Array{Float64,1}) where {T1 <: Number}

Calculates the value of the Husimi Q distribution of a quantum state at the specified points

## args
* rho_in:   d x d array, the density matrix of the state
* xgrid:    vector of floats, x coordinates to calculate the Q function at
* ygrid:    vector of floats, y coordinates to calculate the Q function at

## returns
* Q:        vector of the calculated Q function values
"""
function husimi_Q_pw(rho_in::Array{T1,2},
                     xgrid::Array{Float64,1},
                     ygrid::Array{Float64,1}) where {T1 <: Number}

    size(xgrid,1) == size(ygrid,1) || throw(DimensionMismatch("size of xgrid is not equal to size of ygrid"))

    dim = Int(maximum([3*maximum(abs.(xgrid))^2,3*maximum(abs.(ygrid))^2,size(rho_in,1)]))

    # Disp = zeros(size(rho_in))
    Q = zeros(size(xgrid,1))

    lower, raise, N = create_HO(dim)

    rho_temp = zeros(dim,dim)
    rho_temp = convert(Array{ComplexF64,2},rho_temp)
    rho_temp[1:size(rho_in,1),1:size(rho_in,1)] = rho_in

    state = zeros(dim,1)
    state[1] = 1

    for ii = 1:1:size(xgrid,1)
        @inbounds temp_a = xgrid[ii]+1im*ygrid[ii]
        Disp = exp(temp_a*raise - conj(temp_a)*lower)
        temp_v = Disp*state
        @inbounds Q[ii] = real(temp_v'*rho_temp*temp_v./pi)[1]
    end

    return Q
end

# """
# This function creates the coherent part of a Hamiltonian term in the covariance matrix differential equation.
# Input Arguments:
#     freq: Frequency of the oscillator.
#     sqz: Two-photon (parametric) drive complex amplitude (can be time-dependent)
# """
# function create_h_C(freq::Float64,sqz::Function,t::AbstractFloat)
#     return [imag(sqz(t)) (freq - real(sqz(t)));-(freq + real(sqz(t))) -imag(sqz(t))]
# end


text:
using Dates
using Test
using CliOptions
using CliOptions: consume!

include("testutils.jl")

@testset "Positional()" begin
    @testset "ctor; $(v[1]); $(v[4])" for v in [
        ("single, required", false, nothing, [""], ArgumentError),
        ("single, required", false, nothing, ["-a"], ArgumentError),
        ("single, required", false, nothing, ["a", "-b"], ArgumentError),
        ("single, required", false, nothing, ["a"], ("a",)),
        ("single, required", false, nothing, ["a", "b"], ("a", "b")),
        ("multiple, required", true, nothing, ["a"], ("a",)),
        ("multiple, required", true, nothing, ["a", "b"], ("a", "b")),
        ("single, omittable", false, 42, ["a"], ("a",)),
        ("single, omittable", false, 42, ["a", "b"], ("a", "b")),
    ]
        _, multiple, default, names, expected = v
        if expected isa Type
            @test_throws expected Positional(names...;
                                             multiple = multiple,
                                             default = default)
        else
            option = Positional(names...; multiple = multiple, default = default)
            @test option.names == expected
            @test option.multiple == multiple
            @test option.default === default
        end
    end

    @testset "ctor; duplicates, $(v[1])" for v in [
        (["file", "files"], (nothing, nothing)),
        (["file", "file"], (ArgumentError, "file")),
    ]
        names, expected = v
        if expected[1] isa Type
            tr = @test_throws expected[1] Positional(names...)
            if tr isa Test.Pass
                buf = IOBuffer()
                showerror(buf, tr.value)
                msg = String(take!(buf))
                @test msg == ("ArgumentError: Duplicate names for a Positional found: " *
                              expected[2])
            end
        else
            @test Positional(names...) !== nothing
        end
    end

    @testset "show(x); $(v[1])" for v in [
        (["file"], "Positional(:file)"),
        (["file", "files"], "Positional(:file,:files)"),
    ]
        names, expected_repr = v
        option = Positional(names...)
        @test repr(option) == expected_repr
    end

    @testset "show(io, x)" begin
        let option = Positional("file")
            buf = IOBuffer()
            redirect_stdout(buf) do
                show(option)
            end
            @test String(take!(buf)) == "Positional(:file)"
        end
    end

    @testset "consume!(); $(v[1]), $(v[3])" for v in [
        ("single", false, String[], 0, AssertionError),
        ("single", false, [""], 1, ""),
        ("single", false, ["a"], 1, "a"),
        ("multiple", true, String[], 0, AssertionError),
        ("multiple", true, [""], 1, [""]),
        ("multiple", true, ["a"], 1, ["a"]),
        ("multiple", true, ["a", "b"], 2, ["a", "b"]),
    ]
        _, multiple, args, expected_rv, expected = v
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        option = Positional("file", "files"; multiple = multiple)
        if expected isa Type
            @test_throws expected consume!(d, option, args, ctx)
        else
            next_index = consume!(d, option, args, ctx)
            @test next_index == expected_rv
            @test sorted_keys(d) == ["file", "files"]
            @test d["file"] == expected
        end
    end

    @testset "consume!(); type, $(v[1]), $(v[2])" for v in [
        (Int32, ["2"], Int32(2)),
        (Int32, ["-3"], Int32(-3)),
        (Date, ["2006-01-02"], Date(2006, 1, 2)),
        (Date, ["__not_a_date__"], CliOptionError),
    ]
        T, args, expected = v
        d = Dict{String,Any}()
        ctx = CliOptions.ParseContext()
        option = Positional(T, "value")
        if expected isa Type
            @test_throws expected consume!(d, option, args, ctx)
        else
            consume!(d, option, args, ctx)
            @test d["value"] == expected
        end
    end

    @testset "consume!(); requirement, $(v[1])" for v in [
        ("[7, 13], 13", ["13"],
            Int, [7, 13], (1, 13)),
        ("[7, 13], 7", ["7"],
            Int, [7, 13], (1, 7)),
        ("[7, 13], 0", ["0"],
            Int, [7, 13], (CliOptionError, "must be one of")),
        ("(7, 13), 13", ["13"],
            Int, (7, 13), (1, 13)),
        ("(7, 13), 7", ["7"],
            Int, (7, 13), (1, 7)),
        ("(7, 13), 0", ["0"],
            Int, (7, 13), (CliOptionError, "must be one of")),
        ("/qu+x/, quux", ["quux"],
            String, Regex("qu+x"), (1, "quux")),
        ("/qu+x/, qux", ["qux"],
            String, Regex("qu+x"), (1, "qux")),
        ("/qu+x/, qx", ["qx"],
            String, Regex("qu+x"), (CliOptionError, "must match for")),
        ("String -> Bool, foo", ["foo"],
            String, s -> startswith(s, "foo"), (1, "foo")),
        ("String -> Bool, 6", ["6"],
            Int, n -> iseven(n), (1, 6)),
        ("String -> Bool, 7", ["7"],
            Int, n -> iseven(n), (CliOptionError, "validation failed")),
        ("String -> String, foo", ["foo"],
            String, s -> startswith(s, "foo") ? "" : "It's not foo", (1, "foo")),
        ("String -> String, 6", ["6"],
            Int, n -> iseven(n) ? "" : "must be even", (1, 6)),
        ("String -> String, 7", ["7"],
            Int, n -> iseven(n) ? "" : "must be even", (CliOptionError, "must be even")),
    ]
        _, args, T, requirement, expected = v
        option = Positional(T, "name"; requirement = requirement)
        let d = Dict{String,Any}()
            ctx = CliOptions.ParseContext()
            if expected[1] isa Type
                try
                    CliOptions.consume!(d, option, args, ctx)
                    @test false  # Exception must be thrown
                catch ex
                    @test ex isa expected[1]
                    @test occursin(args[1], ex.msg)
                    @test occursin(expected[2], ex.msg)
                end
            else
                num_consumed = CliOptions.consume!(d, option, args, ctx)
                @test num_consumed == expected[1]
                @test d["name"] == expected[2]
            end
        end
    end

    @testset "check_usage_count(); $(v[1])" for v in [
        ("required, 0", missing, 0, CliOptionError),
        ("required, 1", missing, 1, nothing),
        ("omittable, 0", "foo", 0, nothing),
    ]
        _, default, count, expected = v
        option = Positional("file", default = default)
        ctx = CliOptions.ParseContext()
        ctx.usage_count[option] = count
        if expected isa Type
            @test_throws expected CliOptions.check_usage_count(option, ctx)
        else
            @test true  # No exception was thrown
        end
    end
end


text:
using PGFPlotsX

function μ(x, y, B_p)

    c = B_p.c
    μ_in = B_p.μ_in
    μ_out = B_p.μ_out
    r̄ = B_p.r̄
    r_w = B_p.r_w
    on = B_p.on

    if on == false
        if ndims(x) == 2
            return repeat([μ_out], outer=size(x))
        else
            return repeat([μ_out], outer=length(x))
        end
    else
        return (μ_out - μ_in)/2 *
            (tanh.((x .^ 2 .+ c^2 * y .^ 2 .- r̄) ./ r_w) .+ 1) .+ μ_in
    end
end

function plot_basin(x, y, B_p)

    x = collect(x)
    y = collect(y)
    nx = length(x)
    ny = length(y)

    xm = kron(x, ones(ny)')
    ym = kron(y', ones(nx))

    μ_vals = μ(xm, ym, B_p)

    @pgf basin_plot =  TikzPicture(
    Axis(
        {
            height = "6cm",
            width = "6cm",
            view=(0, 90),
            colorbar_horizontal,
            "colormap/blackwhite",
            colorbar_style={
                xlabel=raw"$\mu$ (GPa)",
                xtick=[B_p.μ_out, B_p.μ_in],
                xticklabels=[raw"$\mu_{out}$", raw"$\mu_{in}$"]
                },
            y_dir="reverse",
            ticks="none",
        },
        Plot3(
            {
                surf,
                shader="flat",
            },
            Coordinates(x,y, μ_vals)
        )
    ),
    Axis(
        {
            height = "6cm",
            width = "6cm",
            axis_lines="box",
            ticks="none",
        }
    ),
    Axis(
        {
            height = "6cm",
            width = "6cm",
            ymin=y[1],
            ymax=y[end],
            xmin=x[1],
            xmax=x[end],
            y_dir="reverse",
            ytick=[4, 8],
            axis_y_line="middle",
            yticklabels=["D", "H"],
            xtick=[-40, 40],
            xticklabels=["-L", "L"],
            xticklabel_pos="top",
        }
        )
    )
    pgfsave("../Basin_paper/figures/basin_plot.tex", basin_plot)

end

function ρ(x, y, B_p)

    c = B_p.c
    ρ_in = B_p.ρ_in
    ρ_out = B_p.ρ_out
    r̄ = B_p.r̄
    r_w = B_p.r_w
    on = B_p.on

    if on == false
        if ndims(x) == 2
            return repeat([ρ_out], outer=size(x))
        else
            return repeat([ρ_out], outer=length(x))
        end
    else
        return (ρ_out - ρ_in)/2 *
            (tanh.((x .^ 2 .+ c^2 * y .^ 2 .- r̄) ./ r_w) .+ 1) .+ ρ_in
    end
end

function μ_x(x, y, B_p)

    c = B_p.c
    μ_in = B_p.μ_in
    μ_out = B_p.μ_out
    r̄ = B_p.r̄
    r_w = B_p.r_w
    on = B_p.on

    if on == false
        if ndims(x) == 2
            return zeros(size(x))
        else
            return zeros(length(x))
        end
    else
        return ((μ_out - μ_in) .* x .*
                sech.((x .^ 2 .+ c^2 * y .^ 2 .- r̄) ./ r_w) .^ 2) ./ r_w
    end
end

function μ_y(x, y, B_p)

    c = B_p.c
    μ_in = B_p.μ_in
    μ_out = B_p.μ_out
    r̄ = B_p.r̄
    r_w = B_p.r_w
    on = B_p.on

    if on == false
        if ndims(x) == 2
            return zeros(size(x))
        else
            return zeros(length(x))
        end

    else
        return ((μ_out - μ_in) .* (c^2 * y) .*
            sech.((x .^ 2 + c^2 * y .^ 2 .- r̄) ./ r_w) .^ 2) ./ r_w
    end
end

function η(y, B_p)
    μf = μ(0, y, B_p)
    return μf ./ (2 .* sqrt.(μf ./ ρ(0, y, B_p)))
end

function fault_params(fc, Dc)

    Wf = 24
    Hvw = 12
    Ht = 6
    δNp = findmin(abs.(Wf .- fc))[2]
    gNp = findmin(abs.(16 .- fc))[2]
    VWp = findmin(abs.((Hvw + Ht) .- fc))[2]
    a = .015
    b0 = .02
    bmin = 0.0

    function b_fun(y)
        if 0 <= y < Hvw
            return b0
        end
        if Hvw <= y < Hvw + Ht
            return b0 + (bmin - b0)*(y-Hvw)/Ht
        end
        if Hvw + Ht <= y < Wf
            return bmin
        end
        return bmin
    end

    RS = (σn = 50.0,
          a = a,
          b = b_fun.(fc),
          Dc = Dc,
          f0 = .6,
          V0 = 1e-6,
          τ_inf = 24.82,
          Vp = 1e-9)

    return δNp, gNp, VWp, RS

end

function plot_ab(a, b, fc)

    a = repeat([a], length(fc))
    @pgf ab_plot = Axis(
        {
            height = "6cm",
            width = "5cm",
            xmin = -0.01,
            xmax = .02,
            xlabel = raw"$a-b$ and $a$",
            ylabel = "Depth (Km)",
            xticklabel_style = ["/pgf/number format/fixed",
                                "/pgf/number format/precision=2",
                                "/pgf/number format/fixed zerofill"],
            scaled_x_ticks = "false",
            legend_entries = [raw"$(a-b)$", raw"$a$"],
            legend_pos = "south west",
            y_dir = "reverse",
        }
    )
    @pgf push!(ab_plot, Plot(
        {
            color = "blue"
        },
        Table(a .- b, fc)
    )
               )
    @pgf push!(ab_plot, Plot(
        {
            color = "red",
            dashed
        },
        Table(a, fc)
    )
               )

    pgfsave("../Basin_paper/figures/ab_plot.tex", ab_plot)

end


text:
using consonance

push!(LOAD_PATH, "/home/ian/Grad/Research/Disks/DiskJockey/")
push!(LOAD_PATH, "/n/home07/iczekala/DiskJockey/")
push!(LOAD_PATH, "/pool/scout0/DiskJockey/")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
end

parsed_args = parse_args(ARGS, s)
run_index = parsed_args["run_index"]

using constants
using image
using visibilities

import PyPlot.plt
using LaTeXStrings

# Log scale the visibility data in intensity
# minv, maxv are passed in true value
function scale(data::Complex128, maxv::Float64)
    if abs(data) > maxv
        return 1.0
    else
        return abs(data)/maxv
    end
end

# Given a 2D matrix, determine the appropriate color scaling
function complex_to_RGB(frame::Matrix{Complex128})
    ny, nx = size(frame)
    HSV = Array(Float64, (ny, nx, 3))
    # min, max = extrema(abs(frame))
    # delta = max - min
    for j=1:ny
        for i=1:nx
            # HSV[j,i,:] = Float64[(angle(frame[j, i]) + pi)/(2pi), 1.0, abs(frame[j, i])/delta]
            HSV[j,i,:] = Float64[(angle(frame[j, i]) + pi)/(2pi), 1.0, scale(frame[j, i], 0.1 * vmax_v)]
        end
    end
    RGB = PyPlot.matplotlib[:colors][:hsv_to_rgb](HSV)
    return RGB
end

global norm = plt.Normalize(vmin=log10(vmax_g) - 6, vmax=log10(vmax_g), clip=false)

function plot_vis(img::image.SkyImage, id::Int)

    fig, ax = plt.subplots(ncols=nchan, nrows=2, figsize=(8, 3.0))

    ext_im = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    for col=1:nchan

        # Plot channel map on top row
        if col != 1
            ax[1,col][:xaxis][:set_ticklabels]([])
            ax[1,col][:yaxis][:set_ticklabels]([])
            ax[1,col][:annotate](@sprintf("%.1f", vels[col]), (0.1, 0.8), xycoords="axes fraction", size=8)
        else
            ax[1,col][:set_xlabel](L"$\Delta \alpha$ ('')", size=8)
            ax[1,col][:set_ylabel](L"$\Delta \delta$ ('')", size=8)
            ax[1,col][:tick_params](axis="both", which="major", labelsize=8)
            ax[1,col][:annotate](@sprintf("%.1f km/s", vels[col]), (0.1, 0.8), xycoords="axes fraction", size=8)
        end

        #Flip the frame for Sky convention
        frame = fliplr(img.data[:,:,col])
        frame += 1e-99 #Add a tiny bit so that we don't have log10(0)
        lframe = log10(frame)
        ax[1,col][:imshow](lframe, extent=ext_im, norm=norm, interpolation="none", origin="lower", cmap=plt.get_cmap("PuBu"))

        # Plot visibilities on bottom

        # Set all labels but the leftmost to blank
        ax[2,col][:tick_params](axis="both", colors="white")
        if col != 1
            ax[2,col][:xaxis][:set_ticklabels]([])
            ax[2,col][:yaxis][:set_ticklabels]([])
            ax[2,col][:xaxis][:set_ticklabels]([])
            ax[2,col][:yaxis][:set_ticklabels]([])
        else
            ax[2,col][:set_xlabel](L"uu [k$\lambda$]", size=8)
            ax[2,col][:set_ylabel](L"vv [k$\lambda$]", size=8)
            ax[2,col][:tick_params](axis="both", which="major", labelsize=8)
            labels = ax[2,col][:get_xticklabels]()
            for label in labels
                label[:set_rotation](40)
                label[:set_color]("black")
            end
            for label in ax[2,col][:get_yticklabels]()
                label[:set_color]("black")
            end
        end

        vis_fft = transform(img, col) # Transform this channel

        # Crop 128 on either side to show more detail
        ncrop = 64
        ext = (vis_fft.uu[end - ncrop], vis_fft.uu[1 + ncrop], vis_fft.vv[1 + ncrop], vis_fft.vv[end - ncrop])
        # Vis needs to be flipped along uu dimension
        frame = fliplr(vis_fft.VV)[1 + ncrop:end-ncrop, 1 + ncrop:end-ncrop]
        ax[2,col][:imshow](complex_to_RGB(frame), interpolation="none", origin="lower", extent=ext)

    end

    # Annotate the parameters
    p = pars[id]
    mass = p.M_star
    r_c = p.r_c
    incl = p.incl
    label = L"$M_\ast$: " * @sprintf("%.2f", mass) * L" $M_\odot$   $r_c$: " * @sprintf("%.0f", r_c) * L" AU   $i$: " * @sprintf("%.0f", incl) * L"${}^\circ$"
    fig[:text](0.35, 0.1, label)

    fig[:subplots_adjust](wspace=0.08, hspace=0.53, top=0.95, bottom=0.26, left=0.1, right=0.9)

    plt.savefig(outdir * @sprintf("vis%04d.png", id), dpi=150)
    plt.close("all")

end

# How many frames per process?
start = run_index * nframes_per_proc + 1

ids = Int[i for i=start:(start + nframes_per_proc)]

for id in ids
    if id <= nframes
        fname = outdir * @sprintf("gimage%04d.out", id)
        im = imread(fname)
        skim = imToSky(im, 73.)
        plot_vis(skim, id)
    else
        break
    end
end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
using Plots, LaTeXStrings, Statistics

function Deposition(;len, tot_time, rate)
    surf = [0 for i=1:len]
    meanList = [0.0 for i=1:tot_time]
    for n in 1:tot_time
        randsurf = rand(1:len,rate)
        for i in randsurf
            surf[i] += 1
        end
        meanList[n] = mean(surf)
    end
    return meanList
end

function Linear_fit(;len, tot_time, rate)
    A = [hcat(1:tot_time) reshape(ones(tot_time), tot_time, 1)]
    b = reshape(meanAvg, tot_time, 1)
    line = (A \ b)
    x = 0:tot_time
    y = x .* line[1] .+ line[2]
    return x, y, line
end

iternum = 1000
Parameters = Dict(
                :len => 200,
                :tot_time => 20,
                :rate => 5000)
allmean = [[0.0 for i in 1:Parameters[:tot_time]] for j = 1:iternum]
meanAvg = [0.0 for i in 1:Parameters[:tot_time]]
means = [0.0 for i in 1:Parameters[:tot_time]]
for i in 1:iternum
    MeanList = Deposition(;Parameters...)
    allmean[i] = MeanList
    meanAvg += MeanList
    print("\r$i")
end
meanAvg /= iternum
for i in 1:Parameters[:tot_time]
    means[i] = std(log.(hcat(allmean...))[i,:])
end

X, Y, Line = Linear_fit(;Parameters...)

scatter(meanAvg,
    xlabel= L"Time",
    ylabel= L"H_{(t)}",
    title= L"H_{(t)}-Time, Particles \ number = 10^{%$(Int(log10(Parameters[:tot_time]*Parameters[:rate])))}",
    yerror = means,
    label = L"Data\ point")
plot!(X,Y,c= :black,label = L"y = %$(round(Line[1],digits= 2))x + %$(round(Line[2],digits= 2))")
savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet2\\Figs\\Q2\\RBDH(t).png")


text:
import Base: convert

"""
An affine map has the general form `y = A*x + b`.

We use `matrix(m)` and `vector(m)` to denote `A` and `b` respectively. Concrete
subtypes include linear maps of the form `y = A*x` and translations of the
form `y = x + b`.
"""
abstract type AbstractAffineMap{T} <: Map{T} end

unsafe_matrix(m::AbstractAffineMap) = m.A
unsafe_vector(m::AbstractAffineMap) = m.b

"Return the matrix `A` in the affine map `Ax+b`."
matrix(m::AbstractAffineMap) = to_matrix(domaintype(m), unsafe_matrix(m), unsafe_vector(m))

"Return the vector `b` in the affine map `Ax+b`."
vector(m::AbstractAffineMap) = to_vector(domaintype(m), unsafe_matrix(m), unsafe_vector(m))

applymap(m::AbstractAffineMap, x) = _applymap(m, x, unsafe_matrix(m), unsafe_vector(m))
_applymap(m::AbstractAffineMap, x, A, b) = A*x + b

applymap!(y, m::AbstractAffineMap, x) = _applymap!(y, m, x, unsafe_matrix(m), unsafe_vector(m))
function _applymap!(y, m::AbstractAffineMap, x, A, b)
    mul!(y, A, x)
    y .+= b
    y
end

isreal(m::AbstractAffineMap) = _isreal(m, unsafe_matrix(m), unsafe_vector(m))
_isreal(m::AbstractAffineMap, A, b) = isreal(A) && isreal(b)

jacobian(m::AbstractAffineMap{T}) where {T} = ConstantMap{T}(matrix(m))
jacobian(m::AbstractAffineMap, x) = matrix(m)

jacdet(m::AbstractAffineMap, x) = _jacdet(m, x, unsafe_matrix(m))
_jacdet(m::AbstractAffineMap, x, A) = det(A)
_jacdet(m::AbstractAffineMap, x::Number, A::UniformScaling) = A.λ
_jacdet(m::AbstractAffineMap, x::AbstractVector, A::Number) = A^length(x)
_jacdet(m::AbstractAffineMap, x::AbstractVector, A::UniformScaling) = A.λ^length(x)

islinear(m::AbstractMap) = false
islinear(m::AbstractAffineMap) = _islinear(m, unsafe_vector(m))
_islinear(m::AbstractAffineMap, b) = all(b .== 0)

isaffine(m::AbstractMap) = islinear(m) || isconstant(m)
isaffine(m::AbstractAffineMap) = true

==(m1::AbstractAffineMap, m2::AbstractAffineMap) =
    matrix(m1) == matrix(m2) && vector(m1) == vector(m2)

==(m1::AbstractAffineMap, m2::IdentityMap) =
    islinear(m1) && matrix(m1) == matrix(m2)
==(m1::IdentityMap, m2::AbstractAffineMap) = m2==m1

mapsize(m::AbstractAffineMap) = _affine_mapsize(m, domaintype(m), unsafe_matrix(m), unsafe_vector(m))
_affine_mapsize(m::AbstractAffineMap, T, A::AbstractArray, b) = size(A)
_affine_mapsize(m::AbstractAffineMap, T, A::AbstractVector, b::AbstractVector) = (length(A),)
_affine_mapsize(m::AbstractAffineMap, T, A::Number, b::Number) = ()
_affine_mapsize(m::AbstractAffineMap, T, A::Number, b::AbstractVector) = (length(b),length(b))
_affine_mapsize(m::AbstractAffineMap, T, A::UniformScaling, b::Number) = ()
_affine_mapsize(m::AbstractAffineMap, T, A::UniformScaling, b) = (length(b),length(b))

Display.displaystencil(m::AbstractAffineMap) = vcat(["x -> "], map_stencil(m, 'x'))
show(io::IO, mime::MIME"text/plain", m::AbstractAffineMap) = composite_show(io, mime, m)

map_stencil(m::AbstractAffineMap, x) = _map_stencil(m, x, unsafe_matrix(m), unsafe_vector(m))
_map_stencil(m::AbstractAffineMap, x, A, b) = [A, " * ", x, " + ", b]
_map_stencil(m::AbstractAffineMap, x, A, b::Real) =
    b >= 0 ? [A, " * ", x, " + ", b] :  [A, " * ", x, " - ", abs(b)]

map_stencil_broadcast(m::AbstractAffineMap, x) = _map_stencil_broadcast(m, x, unsafe_matrix(m), unsafe_vector(m))
_map_stencil_broadcast(m::AbstractAffineMap, x, A, b) = [A, " .* ", x, " .+ ", b]
_map_stencil_broadcast(m::AbstractAffineMap, x, A::Number, b) = [A, " * ", x, " .+ ", b]

Display.object_parentheses(m::AbstractAffineMap) = true
Display.stencil_parentheses(m::AbstractAffineMap) = true

########################
# Linear maps: y = A*x
########################

"""
The supertype of all linear maps `y = A*x`.
Concrete subtypes may differ in how `A` is represented.
"""
abstract type LinearMap{T} <: AbstractAffineMap{T} end

mapsize(m::LinearMap) = _linearmap_size(m, domaintype(m), unsafe_matrix(m))
_linearmap_size(m::LinearMap, T, A) = size(A)
_linearmap_size(m::LinearMap, ::Type{T}, A::Number) where {T<:Number} = ()
_linearmap_size(m::LinearMap, ::Type{T}, A::AbstractVector) where {T<:Number} = (length(A),)
_linearmap_size(m::LinearMap, ::Type{T}, A::Number) where {N,T<:StaticVector{N}} = (N,N)

matrix(m::LinearMap) = to_matrix(domaintype(m), unsafe_matrix(m))
vector(m::LinearMap) = to_vector(domaintype(m), unsafe_matrix(m))

LinearMap(A::Number) = ScalarLinearMap(A)
LinearMap(A::SMatrix) = StaticLinearMap(A)
LinearMap(A::Matrix) = VectorLinearMap(A)
LinearMap(A) = GenericLinearMap(A)

LinearMap{T}(A::Number) where {T<:Number} = ScalarLinearMap{T}(A)
LinearMap{T}(A::SMatrix{M,N}) where {M,N,S,T <: SVector{N,S}} = StaticLinearMap{S}(A)
LinearMap{T}(A::Matrix) where {S,T <: Vector{S}} = VectorLinearMap{S}(A)
LinearMap{T}(A) where {T} = GenericLinearMap{T}(A)

# convenience functions
LinearMap(a::Number...) = LinearMap(promote(a...))
LinearMap(a::NTuple{N,T}) where {N,T <: Number} = LinearMap{SVector{N,T}}(Diagonal(SVector{N,T}(a)))

applymap(m::LinearMap, x) = _applymap(m, x, unsafe_matrix(m))
_applymap(m::LinearMap, x, A) = A*x

applymap!(y, m::LinearMap, x) = _applymap!(y, m, x, unsafe_matrix(m))
_applymap!(y, m::LinearMap, x, A) = mul!(y, A, x)

islinear(m::LinearMap) = true

isreal(m::LinearMap) = _isreal(m, unsafe_matrix(m))
_isreal(m::LinearMap, A) = isreal(A)

==(m1::LinearMap, m2::LinearMap) = matrix(m1) == matrix(m2)

inverse(m::LinearMap) = LinearMap(inv(m.A))
inverse(m::LinearMap, x) = m.A \ x

function leftinverse(m::LinearMap)
    @assert isoverdetermined(m)
    LinearMap(pinv(m.A))
end
function rightinverse(m::LinearMap)
    @assert isunderdetermined(m)
    LinearMap(pinv(m.A))
end
function leftinverse(m::LinearMap, x)
    @assert isoverdetermined(m)
    m.A \ x
end
function rightinverse(m::LinearMap, x)
    @assert isunderdetermined(m)
    m.A \ x
end

similarmap(m::LinearMap, ::Type{T}) where {T} = LinearMap{T}(m.A)

convert(::Type{Map}, a::Number) = LinearMap(a)
convert(::Type{Map{T}}, a::Number) where {T} = LinearMap{T}(a)

convert(::Type{LinearMap}, ::StaticIdentityMap{T}) where {T} = LinearMap{T}(1)
convert(::Type{LinearMap{T}}, ::StaticIdentityMap) where {T} = LinearMap{T}(1)
convert(::Type{AbstractAffineMap}, m::StaticIdentityMap) = convert(LinearMap, m)
convert(::Type{AbstractAffineMap{T}}, m::StaticIdentityMap) where {T} = convert(LinearMap, m)

map_stencil(m::LinearMap, x) = [unsafe_matrix(m), " * ", x]
map_stencil_broadcast(m::LinearMap, x) = _map_stencil_broadcast(m, x, unsafe_matrix(m))
_map_stencil_broadcast(m::LinearMap, x, A) = [A, " .* ", x]
_map_stencil_broadcast(m::LinearMap, x, A::Number) = [A, " * ", x]

"A `GenericLinearMap` is a linear map `y = A*x` for any type of `A`."
struct GenericLinearMap{T,AA} <: LinearMap{T}
    A   ::  AA
end

GenericLinearMap(A) = GenericLinearMap{Any}(A)
GenericLinearMap(A::Number) = GenericLinearMap{typeof(A)}(A)
GenericLinearMap(A::AbstractMatrix{T}) where {T} =
    GenericLinearMap{Vector{T}}(A)
GenericLinearMap(A::SMatrix{M,N,T}) where {M,N,T} =
    GenericLinearMap{SVector{N,T}}(A)
GenericLinearMap(A::AbstractVector{T}) where {T} =
    GenericLinearMap{T}(A)

# Allow any A
GenericLinearMap{T}(A) where {T} = GenericLinearMap{T,typeof(A)}(A)

GenericLinearMap{T}(A::Number) where {T <: Number} = GenericLinearMap{T,T}(A)
GenericLinearMap{T}(A::Number) where {S,T <: AbstractVector{S}} = GenericLinearMap{T,S}(A)
GenericLinearMap{T}(A::AbstractMatrix{S}) where {S,T <: AbstractVector{S}} =
    GenericLinearMap{T,typeof(A)}(A)
GenericLinearMap{T}(A::AbstractMatrix{U}) where {S,T <: AbstractVector{S},U} =
    GenericLinearMap{T}(convert(AbstractMatrix{S}, A))
GenericLinearMap{T}(A::AbstractVector{T}) where {T} =
    GenericLinearMap{T,typeof(A)}(A)
GenericLinearMap{T}(A::AbstractVector{S}) where {S,T} =
    GenericLinearMap{T}(convert(AbstractVector{T}, A))

# Preserve the action on vectors with a number type
inverse(m::GenericLinearMap{T,AA}) where {T<:AbstractVector,AA<:Number} =
    LinearMap{T}(inv(m.A))
leftinverse(m::GenericLinearMap{T,AA}) where {T<:AbstractVector,AA<:Number} =
    LinearMap{T}(inv(m.A))
rightinverse(m::GenericLinearMap{T,AA}) where {T<:AbstractVector,AA<:Number} =
    LinearMap{T}(inv(m.A))

convert(::Type{Map}, A::UniformScaling) = GenericLinearMap{Vector{Any}}(A)
convert(::Type{Map{T}}, A::UniformScaling) where {T} = GenericLinearMap{T}(A)

"A `ScalarLinearMap` is a linear map `y = A*x` for scalars."
struct ScalarLinearMap{T} <: LinearMap{T}
    A   ::  T
end

ScalarLinearMap(A::Number) = ScalarLinearMap{typeof(A)}(A)

isreal(m::ScalarLinearMap{T}) where {T} = isreal(T)

show(io::IO, m::ScalarLinearMap) = show_scalar_linear_map(io, m.A)
show_scalar_linear_map(io, A::Real) = print(io, "x -> $(A) * x")
show_scalar_linear_map(io, A::Complex) = print(io, "x -> ($(A)) * x")

"A `VectorLinearMap` is a linear map `y = A*x` using vectors and matrices."
struct VectorLinearMap{T} <: LinearMap{Vector{T}}
    A   ::  Matrix{T}
end

VectorLinearMap(A::AbstractMatrix{T}) where {T} =
    VectorLinearMap{T}(A)

"A `StaticLinearMap` is a linear map `y = A*x` using static arrays."
struct StaticLinearMap{T,N,M,L} <: LinearMap{SVector{N,T}}
    A   ::  SMatrix{M,N,T,L}
end

StaticLinearMap(A::AbstractMatrix{T}) where {T} =
    StaticLinearMap{T}(A)

StaticLinearMap{T}(A::SMatrix{M,N,S}) where {M,N,T,S} =
    StaticLinearMap{T}(convert(AbstractMatrix{T}, A))
StaticLinearMap{T}(A::SMatrix{M,N,T}) where {M,N,T} =
    StaticLinearMap{T,M,N}(A)
StaticLinearMap{T,N,M}(A::AbstractMatrix) where {T,N,M} =
    StaticLinearMap{T,N,M,M*N}(A)

convert(::Type{Map{SVector{N,T}}}, m::VectorLinearMap) where {N,T} = StaticLinearMap{T,N,N}(m.A)

##########################
# Translations: y = x + b
##########################

"A `Translation` represents the map `y = x + b`."
abstract type Translation{T} <: AbstractAffineMap{T} end

# unsafe_matrix(m::Translation) = I
matrix(m::Translation) = to_matrix(domaintype(m), LinearAlgebra.I, unsafe_vector(m))
vector(m::Translation) = unsafe_vector(m)

mapsize(m::Translation) = _translation_mapsize(m, domaintype(m), unsafe_vector(m))
_translation_mapsize(m::Translation, ::Type{T}, b::Number) where {T} = ()
_translation_mapsize(m::Translation, ::Type{T}, b) where {T} = (length(b),length(b))

map_stencil(m::Translation, x) = _map_stencil(m, x, unsafe_vector(m))
_map_stencil(m::Translation, x, b) = [x, " + ", b]
_map_stencil(m::Translation, x, b::Real) =
    b >= 0 ? [x, " + ", b] : [x, " - ", abs(b)]
map_stencil_broadcast(m::Translation, x) = _map_stencil_broadcast(m, x, unsafe_vector(m))
_map_stencil_broadcast(m::Translation, x, b) = [x, " .+ ", b]

"Translation by a scalar value."
struct ScalarTranslation{T} <: Translation{T}
    b   ::  T
end

isreal(m::ScalarTranslation{T}) where {T} = isreal(T)

show(io::IO, m::ScalarTranslation) = show_scalar_translation(io, m.b)
show_scalar_translation(io, b::Real) = print(io, "x -> x", b < 0 ? " - " : " + ", abs(b))
show_scalar_translation(io, b) = print(io, "x -> x + ", b)

"Translation by a static vector."
struct StaticTranslation{T,N} <: Translation{SVector{N,T}}
    b   ::  SVector{N,T}
end

"Translation by a vector."
struct VectorTranslation{T} <: Translation{Vector{T}}
    b   ::  Vector{T}
end

"Translation by a generic vectorlike object."
struct GenericTranslation{T,B} <: Translation{T}
    b   ::  B
end

Translation(b::Number) = ScalarTranslation(b)
Translation(b::StaticVector) = StaticTranslation(b)
Translation(b::Vector) = VectorTranslation(b)
Translation(b) = GenericTranslation(b)

Translation{T}(b::Number) where {T<:Number} = ScalarTranslation{T}(b)
Translation{T}(b::AbstractVector) where {N,S,T<:SVector{N,S}} = StaticTranslation{S,N}(b)
Translation{T}(b::Vector) where {S,T<:Vector{S}} = VectorTranslation{S}(b)
Translation{T}(b) where {T} = GenericTranslation{T}(b)

jacdet(m::Translation{T}) where {T} = UnityMap{T,eltype(T)}()
jacdet(m::Translation{T}, x) where {T} = one(eltype(T))

isreal(m::Translation) = isreal(unsafe_vector(m))

==(m1::Translation, m2::Translation) = unsafe_vector(m1)==unsafe_vector(m2)

similarmap(m::Translation, ::Type{T}) where {T} = Translation{T}(m.b)

applymap(m::Translation, x) = _applymap(m, x, unsafe_vector(m))
_applymap(m::Translation, x, b) = x + b
applymap!(y, m::Translation, x) = _applymap!(y, m, x, unsafe_vector(m))
_applymap!(y, m::Translation, x, b) = y .= x .+ m.b

inverse(m::Translation{T}) where {T} = Translation{T}(-m.b)
inverse(m::Translation, x) = x - m.b

ScalarTranslation(b::Number) = ScalarTranslation{typeof(b)}(b)

StaticTranslation(b::AbstractVector{T}) where {T} = StaticTranslation{T}(b)

StaticTranslation{T}(b::StaticVector{N}) where {N,T} =
    StaticTranslation{T,N}(b)
StaticTranslation{T}(b::SVector{N,T}) where {N,T} =
    StaticTranslation{T,N}(b)

VectorTranslation(b::AbstractVector{T}) where {T} = VectorTranslation{T}(b)

GenericTranslation(b) = GenericTranslation{typeof(b)}(b)
GenericTranslation(b::AbstractVector{T}) where {T} =
    GenericTranslation{Vector{T}}(b)
GenericTranslation(b::StaticVector{N,T}) where {N,T} =
    GenericTranslation{SVector{N,T}}(b)

GenericTranslation{T}(b) where {T} = GenericTranslation{T,typeof(b)}(b)
GenericTranslation{T}(b::Number) where {T<:Number} =
    GenericTranslation{T,T}(b)

############################
# Affine maps: y = A*x + b
############################

"""
The supertype of all affine maps that store `A` and `b`.
Concrete subtypes differ in how `A` and `b` are represented.
"""
abstract type AffineMap{T} <: AbstractAffineMap{T} end

AffineMap(A::Number, b::Number) = ScalarAffineMap(A, b)
AffineMap(A::SMatrix, b::SVector) = StaticAffineMap(A, b)
AffineMap(A::Matrix, b::Vector) = VectorAffineMap(A, b)
AffineMap(A::UniformScaling{Bool}, b::Number) = ScalarAffineMap(one(b), b)
AffineMap(A, b) = GenericAffineMap(A, b)

AffineMap{T}(A::Number, b::Number) where {T<:Number} = ScalarAffineMap{T}(A, b)
AffineMap{T}(A::AbstractMatrix, b::AbstractVector) where {N,S,T<:SVector{N,S}} = StaticAffineMap{S,N}(A, b)
AffineMap{T}(A::Matrix, b::Vector) where {S,T<:Vector{S}} = VectorAffineMap{S}(A, b)
AffineMap{T}(A::UniformScaling{Bool}, b::Number) where {T} = ScalarAffineMap{T}(one(T), b)
AffineMap{T}(A, b) where {T} = GenericAffineMap{T}(A, b)

similarmap(m::AffineMap, ::Type{T}) where {T} = AffineMap{T}(m.A, m.b)

# If y = A*x+b, then x = inv(A)*(y-b) = inv(A)*y - inv(A)*b
inverse(m::AffineMap) = AffineMap(inv(m.A), -inv(m.A)*m.b)
inverse(m::AffineMap, x) = m.A \ (x-m.b)

# we use matrix_pinv rather than pinv to preserve static matrices
matrix_pinv(A) = pinv(A)
matrix_pinv(A::SMatrix{M,N}) where {M,N} = SMatrix{N,M}(pinv(A))

function leftinverse(m::AffineMap)
    @assert isoverdetermined(m)
    pA = matrix_pinv(m.A)
    AffineMap(pA, -pA*m.b)
end
function rightinverse(m::AffineMap)
    @assert isunderdetermined(m)
    pA = matrix_pinv(m.A)
    AffineMap(pA, -pA*m.b)
end
function leftinverse(m::AffineMap, x)
    @assert isoverdetermined(m)
    m.A \ (x-m.b)
end
function rightinverse(m::AffineMap, x)
    @assert isunderdetermined(m)
    m.A \ (x-m.b)
end

"An affine map for any combination of types of `A` and `b`."
struct GenericAffineMap{T,AA,B} <: AffineMap{T}
    A   ::  AA
    b   ::  B
end

GenericAffineMap(A, b) = GenericAffineMap{typeof(b)}(A, b)
GenericAffineMap(A::AbstractVector{S}, b::AbstractVector{T}) where {S,T} =
    GenericAffineMap{promote_type(S,T)}(A, b)
GenericAffineMap(A::AbstractArray{S}, b::AbstractVector{T}) where {S,T} =
    GenericAffineMap{Vector{promote_type(S,T)}}(A, b)
GenericAffineMap(A::SMatrix{M,N,S}, b::StaticVector{M,T}) where {M,N,S,T} =
    GenericAffineMap{SVector{N,promote_type(S,T)}}(A, b)
GenericAffineMap(A::SMatrix{M,N,S}, b::AbstractVector{T}) where {M,N,S,T} =
    GenericAffineMap{SVector{N,promote_type(S,T)}}(A, b)
GenericAffineMap(A::S, b::AbstractVector{T}) where {S<:Number,T} =
    GenericAffineMap{Vector{promote_type(S,T)}}(A, b)
GenericAffineMap(A::S, b::StaticVector{N,T}) where {S<:Number,N,T} =
    GenericAffineMap{SVector{N,promote_type(S,T)}}(A, b)
GenericAffineMap(A::UniformScaling{Bool}, b) =
    GenericAffineMap(UniformScaling{eltype(b)}(1), b)

# Fallback routine for generic A and b, special cases follow
GenericAffineMap{T}(A, b) where {T} = GenericAffineMap{T,typeof(A),typeof(b)}(A, b)

GenericAffineMap{T}(A::AbstractVector{S}, b::AbstractVector{U}) where {T<:Number,S,U} =
    GenericAffineMap{T}(convert(AbstractVector{T}, A), convert(AbstractVector{T}, b))
GenericAffineMap{T}(A::AbstractVector{T}, b::AbstractVector{T}) where {T<:Number} =
    GenericAffineMap{T,typeof(A),typeof(b)}(A, b)
GenericAffineMap{T}(A::Number, b) where {T} = GenericAffineMap{T,eltype(T),typeof(b)}(A, b)
GenericAffineMap{T}(A::Number, b::AbstractVector) where {N,S,T <: SVector{N,S}} =
    GenericAffineMap{T,S,SVector{N,S}}(convert(S,A), SVector{N,S}(b))
# Promote element types of abstract arrays
GenericAffineMap{T}(A::AbstractArray, b::AbstractVector) where {S,T<:AbstractVector{S}} =
    GenericAffineMap{T}(convert(AbstractArray{eltype(T)},A), convert(AbstractVector{eltype(T)}, b))
GenericAffineMap{T}(A::AbstractArray{S}, b::AbstractVector{S}) where {S,T<:AbstractVector{S}} =
    GenericAffineMap{T,typeof(A),typeof(b)}(A, b)
GenericAffineMap{T}(A::UniformScaling{Bool}, b::AbstractVector) where {S,T<:AbstractVector{S}} =
    GenericAffineMap{T}(A*one(S), convert(AbstractVector{S}, b))
GenericAffineMap{T}(A::UniformScaling{S}, b::AbstractVector{S}) where {S,T<:AbstractVector{S}} =
    GenericAffineMap{T,typeof(A),typeof(b)}(A, b)

similarmap(m::GenericAffineMap, ::Type{T}) where {T} = AffineMap{T}(m.A, m.b)

convert(::Type{GenericAffineMap{T}}, m::GenericAffineMap) where {T} =
    GenericAffineMap{T}(m.A, m.b)

"An affine map with scalar representation."
struct ScalarAffineMap{T} <: AffineMap{T}
    A   ::  T
    b   ::  T
end

ScalarAffineMap(A, b) = ScalarAffineMap(promote(A, b)...)

isreal(m::ScalarAffineMap{T}) where {T} = isreal(T)

show(io::IO, m::ScalarAffineMap) = show_scalar_affine_map(io, m.A, m.b)
show_scalar_affine_map(io, A::Real, b::Real) = print(io, "x -> $(A) * x", b < 0 ? " - " : " + ", abs(b))
show_scalar_affine_map(io, A::Complex, b::Complex) = print(io, "x -> ($(A)) * x + ", b)
show_scalar_affine_map(io, A, b) = print(io, "x -> ($(A)) * x + $(b)")

convert(::Type{ScalarAffineMap{T}}, m::ScalarAffineMap) where {T} =
    ScalarAffineMap{T}(m.A, m.b)

"An affine map with array and vector representation."
struct VectorAffineMap{T} <: AffineMap{Vector{T}}
    A   ::  Matrix{T}
    b   ::  Vector{T}
end

VectorAffineMap(A::AbstractArray{T}, b::AbstractVector{T}) where {T} =
    VectorAffineMap{T}(A, b)
function VectorAffineMap(A::AbstractArray{S}, b::AbstractVector{T}) where {S,T}
    U = promote_type(S,T)
    VectorAffineMap(convert(AbstractArray{U}, A), convert(AbstractVector{U}, b))
end

convert(::Type{VectorAffineMap{T}}, m::VectorAffineMap) where {T} =
    VectorAffineMap{T}(m.A, m.b)

"An affine map with representation using static arrays."
struct StaticAffineMap{T,N,M,L} <: AffineMap{SVector{N,T}}
    A   ::  SMatrix{M,N,T,L}
    b   ::  SVector{M,T}
end

# Constructors:
# - first, we deduce T
StaticAffineMap(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T} =
    StaticAffineMap{T}(A, b)
function StaticAffineMap(A::AbstractMatrix{S}, b::AbstractVector{T}) where {S,T}
    U = promote_type(S,T)
    StaticAffineMap(convert(AbstractMatrix{U}, A), convert(AbstractVector{U}, b))
end

StaticAffineMap{T}(A::AbstractMatrix, b::AbstractVector) where {T} =
    StaticAffineMap{T}(convert(AbstractMatrix{T}, A), convert(AbstractVector{T}, b))

# - then, we determine N and/or M, from the arguments
function StaticAffineMap{T}(A::AbstractMatrix{T}, b::SVector{M,T}) where {T,M}
    @assert size(A) == (M,M)
    StaticAffineMap{T,M,M}(A, b)
end
StaticAffineMap{T}(A::SMatrix{M,N,T}, b::AbstractVector) where {T,N,M} =
    StaticAffineMap{T,N,M}(A, b)
StaticAffineMap{T}(A::SMatrix{M,N,T}, b::SVector{M,T}) where {T,N,M} =
    StaticAffineMap{T,N,M}(A, b)
StaticAffineMap{T,N}(A::AbstractMatrix, b::AbstractVector) where {T,N} =
    StaticAffineMap{T,N,N}(A, b)
StaticAffineMap{T,N}(A::SMatrix{M,N}, b::AbstractVector) where {T,N,M} =
    StaticAffineMap{T,N,M}(A, b)

# - finally invoke the constructor (and implicitly convert the data if necessary)
StaticAffineMap{T,N,M}(A::AbstractMatrix, b::AbstractVector) where {T,N,M} =
    StaticAffineMap{T,N,M,M*N}(A, b)

convert(::Type{Map{SVector{N,T}}}, m::VectorAffineMap) where {N,T} =
    StaticAffineMap{T,N}(m.A, m.b)


text:
# subroutine to compute a default number of folds
@inline cv_get_num_folds(nmin::Int, nmax::Int) = max(nmin, min(Sys.CPU_CORES::Int, nmax))

"A subroutine to efficiently compute i = indmax(abs(x))."
function indmax_abs{T <: Float}(x::DenseVector{T})
    a = abs(x[1])
    b = 1
    for i in eachindex(x)
        c = abs(x[i])
        if c > a
            a = c
            b = i
        end
    end
    return b
end

# verbose printing of cv results
function print_cv_results{T <: Float}(io::IO, errors::Vector{T}, path::DenseVector{Int}, k::Int)
    println(io, "\n\nCrossvalidation Results:")
    println(io, "k\tMSE")
    for i = 1:length(errors)
        println(io, path[i], "\t", errors[i])
    end
    println(io, "\nThe lowest MSE is achieved at k = ", k)
end

# default IO for print_cv_results is STDOUT
print_cv_results{T <: Float}(errors::Vector{T}, path::DenseVector{Int}, k::Int) = print_cv_results(STDOUT, errors, path, k)


text:
type Stream
  ptr::cudaStream_t
end

function Stream()
  p = cudaStream_t[0]
  cudaStreamCreate(p)
  p = p[1]
  Stream(p)
end

const nullstream = Stream()

Base.unsafe_convert(::Type{cudaStream_t}, s::Stream) = s.ptr

destroy(s::Stream) = cudaStreamDestroy(s)

synchronize(s::Stream) = cudaStreamSynchronize(s)