load_file("Cell_positions.hoc")
load_file("Rod.tem")
load_file("Cone.tem")
load_file("Horizontal.tem")
load_file("Bip.tem")
load_file("A2.tem")
load_file("nrngui.hoc")
load_file("read_population.hoc")
load_file("Parameters.hoc")

calc_cell_density(0, 0, 0, 1)
calc_cell_spacing()
determine_cell_init_positions()
generate_population_positions()

PR_total_cells = 0
    rod_total_cells = 0
    cone_total_cells = 0
HZ_total_cells = 0
BIP_total_cells = 0
    rod_bip_total_cells = 0
    ON_bip_total_cells = 0
    OFF_bip_total_cells = 0
AII_total_cells = 0

PR_total_upper_cells = 0
PR_total_lower_cells = 0
HZ_total_upper_cells = 0
HZ_total_lower_cells = 0
BIP_total_upper_cells = 0
BIP_total_lower_cells = 0
AII_total_upper_cells = 0
AII_total_lower_cells = 0

// objects for cells
objref PR_upper[PR_count_upper_cells], PR_lower[PR_count_lower_cells - PR_total_diag_cells]
objref HZ_upper[HZ_count_upper_cells], HZ_lower[HZ_count_lower_cells - HZ_total_diag_cells]
objref BIP_upper[BIP_count_upper_cells], BIP_lower[BIP_count_lower_cells - BIP_total_diag_cells]
objref AII_upper[AII_count_upper_cells], AII_lower[AII_count_lower_cells - AII_total_diag_cells]

// objects for saving cell positions
objref PR_upper_positions, PR_lower_positions
objref HZ_upper_positions, HZ_lower_positions
objref BIP_upper_positions, BIP_lower_positions
objref AII_upper_positions, AII_lower_positions

// objects for saving cell index
objref PR_upper_index, PR_lower_index
objref HZ_upper_index, HZ_lower_index
objref BIP_upper_index, BIP_lower_index
objref AII_upper_index, AII_lower_index

// objects for saving type
objref PR_upper_cell_type, PR_lower_cell_type
objref BIP_upper_cell_type, BIP_lower_cell_type

objref line_for_display

// initialising matrices for saving positions
PR_upper_positions = new Matrix(PR_count_upper_cells, 2)
PR_lower_positions = new Matrix(PR_count_lower_cells - PR_total_diag_cells, 2)
HZ_upper_positions = new Matrix(HZ_count_upper_cells, 2)
HZ_lower_positions = new Matrix(HZ_count_lower_cells - HZ_total_diag_cells, 2)
BIP_upper_positions = new Matrix(BIP_count_upper_cells, 2)
BIP_lower_positions = new Matrix(BIP_count_lower_cells - BIP_total_diag_cells, 2)
AII_upper_positions = new Matrix(AII_count_upper_cells, 2)
AII_lower_positions = new Matrix(AII_count_lower_cells - AII_total_diag_cells, 2)

// initialising matrices for saving cell indices
// First value is the index position, Second value is the subtype of the cell
PR_upper_index = new Matrix(PR_count_upper_cells, 1)
PR_lower_index = new Matrix(PR_count_lower_cells - PR_total_diag_cells, 1)
HZ_upper_index = new Matrix(HZ_count_upper_cells, 1)
HZ_lower_index = new Matrix(HZ_count_lower_cells - HZ_total_diag_cells, 1)
BIP_upper_index = new Matrix(BIP_count_upper_cells, 1)
BIP_lower_index = new Matrix(BIP_count_lower_cells - BIP_total_diag_cells, 1)
AII_upper_index = new Matrix(AII_count_upper_cells, 1)
AII_lower_index = new Matrix(AII_count_lower_cells - AII_total_diag_cells, 1)

PR_upper_cell_type = new Matrix(PR_count_upper_cells, 1)
PR_lower_cell_type = new Matrix(PR_count_lower_cells - PR_total_diag_cells, 1)
BIP_upper_cell_type = new Matrix(BIP_count_upper_cells, 1)
BIP_lower_cell_type = new Matrix(BIP_count_lower_cells - BIP_total_diag_cells, 1)

proc generate_PR_array() {
    // Generate the diagonal cells and upper half
    k = 0
    for i = 0, PR_total_diag_cells - 1 {
        // Determine whether cell is rod or cone based on distribution
        for j = 0, PR_count_upper.x[i] - 1 {
            calc_PR_density(PR_population_upper_position[i].x[j][0], PR_population_upper_position[i].x[j][1])
            if(abs(PR_population_upper_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                if(PR_type_value == 1) {
                    PR_upper[k] = new Rod(PR_population_upper_position[i].x[j][0]*1000, PR_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness + OPL_Thickness + ONL_PIS_Thickness)*1000)
                    PR_upper_positions.x[k][0] = PR_population_upper_position[i].x[j][0]*1000
                    PR_upper_positions.x[k][1] = PR_population_upper_position[i].x[j][1]*1000
                    PR_upper_index.x[k][0] = i + j
                    PR_upper_cell_type.x[k][0] = PR_type_value
                    k = k + 1
                    rod_total_cells = rod_total_cells + 1
                }
                if(PR_type_value == 2) {
                    PR_upper[k] = new Cone(PR_population_upper_position[i].x[j][0]*1000, PR_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness + OPL_Thickness + ONL_PIS_Thickness)*1000)
                    PR_upper_positions.x[k][0] = PR_population_upper_position[i].x[j][0]*1000
                    PR_upper_positions.x[k][1] = PR_population_upper_position[i].x[j][1]*1000
                    PR_upper_index.x[k][0] = i + j
                    PR_upper_cell_type.x[k][0] = PR_type_value
                    k = k + 1
                    cone_total_cells = cone_total_cells + 1
                }
            }
        }
    }
    PR_total_upper_cells = k

    // Generate the lower half
    l = 0
    for i = 0, PR_total_diag_cells - 1 {
        // we use j = 1 to skip the diagonal cells already generated by the upper half
        for j = 1, PR_count_lower.x[i] - 1 {
            calc_PR_density(PR_population_lower_position[i].x[j][0], PR_population_lower_position[i].x[j][1])
            if(abs(PR_population_lower_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                if(PR_type_value == 1) {
                    PR_lower[l] = new Rod(PR_population_lower_position[i].x[j][0]*1000, PR_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness + OPL_Thickness + ONL_PIS_Thickness)*1000)
                    PR_lower_positions.x[l][0] = PR_population_lower_position[i].x[j][0]*1000
                    PR_lower_positions.x[l][1] = PR_population_lower_position[i].x[j][1]*1000
                    PR_lower_index.x[l][0] = i - j
                    PR_lower_cell_type.x[l][0] = PR_type_value
                    l = l + 1
                    rod_total_cells = rod_total_cells + 1
                }
                if(PR_type_value == 2) {
                    PR_lower[l] = new Cone(PR_population_lower_position[i].x[j][0]*1000, PR_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness + OPL_Thickness + ONL_PIS_Thickness)*1000)
                    PR_lower_positions.x[l][0] = PR_population_lower_position[i].x[j][0]*1000
                    PR_lower_positions.x[l][1] = PR_population_lower_position[i].x[j][1]*1000
                    PR_lower_index.x[l][0] = i - j
                    PR_lower_cell_type.x[l][0] = PR_type_value
                    l = l + 1
                    cone_total_cells = cone_total_cells + 1
                }
            }
        }
    } 
    PR_total_lower_cells = l

    PR_total_cells = rod_total_cells + cone_total_cells
}

proc generate_HZ_array() {
    // Generate the diagonal cells and upper half
    k = 0
    for i = 0, HZ_total_diag_cells - 1 {
        // Determine whether cell is rod or cone based on distribution
        for j = 0, HZ_count_upper.x[i] - 1 {
            if(abs(HZ_population_upper_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                if(abs(HZ_population_upper_position[i].x[j][0] - Eccentricity_x) < (Patch_size_x/2)) {
                    HZ_upper[k] = new Horizontal(HZ_population_upper_position[i].x[j][0]*1000, HZ_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness + OPL_Thickness)*1000)
                    HZ_upper_positions.x[k][0] = HZ_population_upper_position[i].x[j][0]*1000
                    HZ_upper_positions.x[k][1] = HZ_population_upper_position[i].x[j][1]*1000
                    HZ_upper_index.x[k][0] = i + j
                    k = k + 1
                    HZ_total_cells = HZ_total_cells + 1
                }
            }
        }
    }
    HZ_total_upper_cells = k

    // Generate the lower half
    l = 0
    for i = 0, HZ_total_diag_cells - 1 {
        // we use j = 1 to skip the diagonal cells already generated by the upper half
        for j = 1, HZ_count_lower.x[i] - 1 {
            if(abs(HZ_population_lower_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                if(abs(HZ_population_lower_position[i].x[j][0] - Eccentricity_x) < (Patch_size_x/2)) {
                    HZ_lower[l] = new Horizontal(HZ_population_lower_position[i].x[j][0]*1000, HZ_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness + OPL_Thickness)*1000)
                    HZ_lower_positions.x[l][0] = HZ_population_lower_position[i].x[j][0]*1000
                    HZ_lower_positions.x[l][1] = HZ_population_lower_position[i].x[j][1]*1000
                    HZ_lower_index.x[l][0] = i - j
                    l = l + 1
                    HZ_total_cells = HZ_total_cells + 1
                }
            }
        }
    }
    HZ_total_lower_cells = l
    
}

proc generate_BIP_array() {
    // Generate the diagonal cells and upper half
    k = 0
    for i = 0, BIP_total_diag_cells - 1 {
        // Determine whether cell is RBC, ON CBC or OFF CBC based on distribution
        for j = 0, BIP_count_upper.x[i] - 1 {
            calc_BIP_density(BIP_population_upper_position[i].x[j][0], PR_population_upper_position[i].x[j][1])
            if(abs(BIP_population_upper_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                if(BIP_type_value == 1) {
                    BIP_upper[k] = new Bip(BIP_population_upper_position[i].x[j][0]*1000, BIP_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness)*1000)
                    BIP_upper_positions.x[k][0] = BIP_population_upper_position[i].x[j][0]*1000
                    BIP_upper_positions.x[k][1] = BIP_population_upper_position[i].x[j][1]*1000
                    BIP_upper_index.x[k][0] = i + j
                    BIP_upper_cell_type.x[k][0] = BIP_type_value
                    k = k + 1
                    rod_bip_total_cells = rod_bip_total_cells + 1
                }
                if(BIP_type_value == 2) {
                    BIP_upper[k] = new Bip(BIP_population_upper_position[i].x[j][0]*1000, BIP_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness)*1000)
                    BIP_upper_positions.x[k][0] = BIP_population_upper_position[i].x[j][0]*1000
                    BIP_upper_positions.x[k][1] = BIP_population_upper_position[i].x[j][1]*1000
                    BIP_upper_index.x[k][0] = i + j
                    BIP_upper_cell_type.x[k][0] = BIP_type_value
                    k = k + 1
                    ON_bip_total_cells = ON_bip_total_cells + 1
                }
                if(BIP_type_value == 3) {
                    BIP_upper[k] = new Bip(BIP_population_upper_position[i].x[j][0]*1000, BIP_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness)*1000)
                    BIP_upper_positions.x[k][0] = BIP_population_upper_position[i].x[j][0]*1000
                    BIP_upper_positions.x[k][1] = BIP_population_upper_position[i].x[j][1]*1000
                    BIP_upper_index.x[k][0] = i + j
                    BIP_upper_cell_type.x[k][0] = BIP_type_value
                    k = k + 1
                    OFF_bip_total_cells = OFF_bip_total_cells + 1
                }
            }
        }
    }
    BIP_total_upper_cells = k

    // Generate the lower half
    l = 0
    for i = 0, BIP_total_diag_cells - 1 {
        // we use j = 1 to skip the diagonal cells already generated by the upper half
        for j = 1, BIP_count_lower.x[i] - 1 {
            calc_BIP_density(BIP_population_lower_position[i].x[j][0], BIP_population_lower_position[i].x[j][1])
            if(abs(BIP_population_lower_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                if(BIP_type_value == 1) {
                    BIP_lower[l] = new Bip(BIP_population_lower_position[i].x[j][0]*1000, BIP_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness)*1000)
                    BIP_lower_positions.x[l][0] = BIP_population_lower_position[i].x[j][0]*1000
                    BIP_lower_positions.x[l][1] = BIP_population_lower_position[i].x[j][1]*1000
                    BIP_lower_index.x[l][0] = i - j
                    BIP_lower_cell_type.x[l][0] = BIP_type_value
                    l = l + 1
                    rod_bip_total_cells = rod_bip_total_cells + 1
                }
                if(BIP_type_value == 2) {
                    BIP_lower[l] = new Bip(BIP_population_lower_position[i].x[j][0]*1000, BIP_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness)*1000)
                    BIP_lower_positions.x[l][0] = BIP_population_lower_position[i].x[j][0]*1000
                    BIP_lower_positions.x[l][1] = BIP_population_lower_position[i].x[j][1]*1000
                    BIP_lower_index.x[l][0] = i - j
                    BIP_lower_cell_type.x[l][0] = BIP_type_value
                    l = l + 1
                    ON_bip_total_cells = ON_bip_total_cells + 1
                }
                if(BIP_type_value == 3) {
                    BIP_lower[l] = new Bip(BIP_population_lower_position[i].x[j][0]*1000, BIP_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness + INL_Thickness)*1000)
                    BIP_lower_positions.x[l][0] = BIP_population_lower_position[i].x[j][0]*1000
                    BIP_lower_positions.x[l][1] = BIP_population_lower_position[i].x[j][1]*1000
                    BIP_lower_index.x[l][0] = i - j
                    BIP_lower_cell_type.x[l][0] = BIP_type_value
                    l = l + 1
                    OFF_bip_total_cells = OFF_bip_total_cells + 1
                }
            }
        }
    }
    BIP_total_lower_cells = l

    BIP_total_cells = rod_bip_total_cells + ON_bip_total_cells + OFF_bip_total_cells
}

proc generate_AII_array() {
    // Generate the diagonal cells and upper half
    k = 0
    for i = 0, AII_total_diag_cells - 1 {
        // Determine whether cell is rod or cone based on distribution
        for j = 0, AII_count_upper.x[i] - 1 {
            if(abs(AII_population_upper_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                AII_upper[k] = new A2(AII_population_upper_position[i].x[j][0]*1000, AII_population_upper_position[i].x[j][1]*1000, (IPL_GCL_Thickness)*1000)
                AII_upper_positions.x[k][0] = AII_population_upper_position[i].x[j][0]*1000
                AII_upper_positions.x[k][1] = AII_population_upper_position[i].x[j][1]*1000
                AII_upper_index.x[k][0] = i + j
                k = k + 1
                AII_total_cells = AII_total_cells + 1
            }
        }
    }
    AII_total_upper_cells = k

    // Generate the lower half
    l = 0
    for i = 0, AII_total_diag_cells - 1 {
        // we use j = 1 to skip the diagonal cells already generated by the upper half
        for j = 1, AII_count_lower.x[i] - 1 {
            if(abs(AII_population_lower_position[i].x[j][1] - Eccentricity_y) < (Patch_size_y/2)) {
                AII_lower[l] = new A2(AII_population_lower_position[i].x[j][0]*1000, AII_population_lower_position[i].x[j][1]*1000, (IPL_GCL_Thickness)*1000)
                AII_lower_positions.x[l][0] = AII_population_lower_position[i].x[j][0]*1000
                AII_lower_positions.x[l][1] = AII_population_lower_position[i].x[j][1]*1000
                AII_lower_index.x[l][0] = i - j
                l = l + 1
                AII_total_cells = AII_total_cells + 1
            }
        }
    } 
    AII_total_lower_cells = l

}

proc generate_RGC_array() {

    read_write_RGC_population()
    split_RGC_parameters()
    load_RGC_templates()
    generate_RGC_cells()
    calc_dendrite_centre()
    change_RGC_ionic_properties()

    //line_for_display = new Line(-1200, 0, 100)

}

proc print_total_cells() {
    // Procedure for printing the total number of cells
    print "Total number of PRs is"
    print PR_total_cells
    print "Total number of rods is"
    print rod_total_cells
    print "Total number of cones is"
    print cone_total_cells
    print "Total number of HZs is"
    print HZ_total_cells
    print "Total number of BIPs is"
    print BIP_total_cells
    print "Total number of rod bips is"
    print rod_bip_total_cells
    print "Total number of ON bips is"
    print ON_bip_total_cells
    print "Total number of OFF bips is"
    print OFF_bip_total_cells
    print "Total number of AII ACs is"
    print AII_total_cells
}